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ABSTRACT 


A least squares algorithm is developed to solve for the 
trajectory and transponder array coordinates of the current 
velocity profiler, Pegasus. Measurement residuals and 
parameter precision are computed for data quality analysis. 
Travel times from a maximum of four seafloor transponders, 
pressure sensor depths, and transponder positions are input 
with their respective accuracy estimates. The algorithm is 
used to analyze a 2250 m profile from the Monterey Canyon with 
four transponders, one of which had not been positioned. This 
transponder's unknown position is found and problems in the 
other array coordinates identified. Transponder coordinate 
precision improves by factors of ten in the horizontal and 
five in depth, to about 13 m (Drms) and 2 m (10) respectively. 
Trajectory precision is about 7 m (Drms) horizontal, with high 
correlation between points. Thus, the precision of horizontal 
velocity components, determined by time differencing points, 
is better than 11 cm/s (10). Depth precision is better than 
3 m (lo), except in the deeper portions where anomalous 
pressure residuals near the depth of the transponder array 


suggest systematic pressure errors needing further study. 
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I. INTRODUCTION 


The Oceanography Department of the Naval Postgraduate 
School is presently conducting a study of ocean current 
velocities in the waters off the central California coast in 
the area of the Monterey Canyon. During the course of the 
study, several stations have been established for the 
collection of data. Station 10, northernmost in this chain of 
stations and the primary focus for analyses presented in this 
thesis, contains four seafloor transponders, one of which was 
considered to have an essentially unknown position at the 
beginning of the work described here. 

The main oceanographic instrument used in this survey is 
the dropsonde Pegasus. As it descends and ascends through the 
water column, it records at 16 second intervals, both the data 
from oceanographic sensors and the time for return travel of 
acoustic signals from transponders previously set in place on 
the ocean floor. 

In the past, Pegasus-generated velocity studies have used 
the raw data in various ways. The trajectory of the 
instrument has been fitted to a curve which was then differen- 
tiated with respect to time in order to produce the horizontal 
velocity components (Halkin et al., 1985). Alternatively, and 


as is the present practice at NPS, velocities have been 


calculated from position differences and then smoothed with a 
spatial filter such as a seven point running average. 

This study outlines a method of determining the current 
velocities and their precision by applying a least squares 
adjustment procedure, adapted from geodetic survey techniques, 
whereby all observational data are used simultaneously. 
Specifically, three basic problems will be analyzed: 


* The precision of the transponder coordinates from which 
the position of Pegasus is derived. 


x The precision of the Pegasus positions and consequently 
the current velocities derived from Pegasus navigation. 


* The depth at which pressure becomes critical in solving 
for velocities. 


Discussions of this topic will appear in the following 
sequence: 

* Chapter II provides background information on the 

transponder network and the Pegasus system and its 


positioning. 


* Chapter III discusses sources of error resulting from 
systematic and observational inaccuracies in the data. 


* Chapter IV gives a brief explanation of the least squares 
adjustment process. 


* Chapter V describes the final procedure used. 
* Chapter VI discusses the results obtained. 


* Chapter VII provides conclusions and recommendations for 
future consideration. 


* Appendices contain algorithms and explanations of the 
Fortran programs used. 


II. BACKGROUND 


The area being surveyed iS composed of ten stations 
located just off the California coast approximately between 
36°05’ and 36°39’ North Latitude, and 122°08’ and 124°13’ West 
Longitude. Of particular interest here is Station C10, 
northernmost in the chain and located in the Monterey Canyon 


about ten miles west of Monterey (see Figure 1). 


A. TRANSPONDER NETWORK 

At this station four transponders were Aenea, horizon- 
tally separated by distances somewhat equivalent to transpon- 
der depths (around 2000 m). Three of these transponders, 
responding at 11.5 kHz, 12.0 kHz, and 12.5 kHz, respectively, 
were located by the survey vessel according to traditional 
methods as described below. However, the bandwidth of the 
shipboard receiving instrument was too narrow to pick up the 
13.5 kHz signal from the fourth transponder, and therefore its 
position was largely unknown. (Table 1 defines transponder 
abbreviations, their frequencies, and their approximate 


depths. Figure 2 shows the survey net for Station C10.) 
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Figure 1. General Area of Station C10 


TABLE 1 
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Figure 2. Station C10 Survey Net 
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B. SURVEYING A TRANSPONDER NETWORK 
1. Transponder Depth 
To determine each available transponder depth, the 
Survey ship homes in on the instrument's signal until a 
minimum time 1s recorded, thereby assuming the ship to be 
directly over the instrument (see Figure 3). An appropriate 
harmonic mean sound velocity multiplied by one way signal 


travel time produces the depth of the transponder. 





SAO 
TRACK 
ii 
‘jj TRANSPONDER 
) 
Figure 3. Determining Transponder Depth 
2. Baselines 


To establish the relative horizontal x and y coordi- 
nates of transponders, baselines are determined between the 
instruments. Running at a relatively slow speed (say 4-6 


kts), the survey vessel attempts to cross each baseline at a 


90° angle. The process is repeated in the opposite direction. 
This procedure is then duplicated, several more times if 
possible. Gareful ‘observation of the analog trace that 
records round-trip travel times of the acoustic signal between 
the survey vessel and the two transponders indicates when the 


shortest distance has been reached, as shown in Figure 4. 
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Figure 4. Analog Trace of a Baseline Crossing 


Signals show up on the trace as nena leone approaching 
one another. They will lie vertically in line if the baseline 
is crossed at a 90° angle and therefore the shortest distance to 
each of them is reached Simultaneously. At this point the 
ship will lie in a vertical plane containing the two 


transponders. The sum of the two horizontal ranges will be 


at a minimum and therefore establish the baseline (see Figure 





5). 
SHIP TRACK 
Figure 5. Baseline Determination 
The time of crossing, ship's course, latitude/ 
longitude, and Loran coordinates are noted. Time units for 


each minimum range are measured and then multiplied by the 
appropriate harmonic mean sound velocity to obtain slope 
distances. These distances, the depths of the two transpon- 
ders, and the ship's course are then used to calculate the 
horizontal baseline length and the azimuth between the two 
transponders. Finally, a local xy plane coordinate system is 
established, setting the position of one of the transponders 
at 0,0, and positioning the other transponder relative to it. 

Many factors combine to make the information on this 


analog trace subjective and inaccurate. The analog paper may 


be set to progress at different rates producing sweeps of 
eos) | Sewer 2 Ssyeete. The trace may actually "wrap around" 
one or more times as the time interval gets longer and longer, 
making page integer resolution difficult. At a one second 
Sweep, the full page represents (after conversion from time to 
distance) meters traveled in one second, i.e., 1500 m/s. 
Furthermore, if the baseline is not crossed perpendic- 
ularly, a time offset between minimum ranges is observed which 
decreases the length of the baseline. This error can be 
removed to some extent by multiplying the time offset by the 
ship's speed and then by geometric relationships, computing 
ene Correction. This situation is illustrated in Figure 6 
where T1-T2 represents the actual baseline, A, B represent 
respective points on the analog trace where minimum slant 
ranges are indicated, and a, b represent respective minimum 


slant ranges to each transponder. 
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Figure 6. Time Offset 


In addition to the above possibilities for error, an 
inaccurate course heading at the time of baseline crossing 
will alter the geometry of the array, resulting in inaccurate 
coordinates. Further discussions of transponder depth and 
baseline determinations can be found in Kuo (1985), McKeown 
(1975), and Hart V1967)< For an assessment of transponder 
network accuracy, refer to Chapter III. 

Figure 7 illustrates the C10 transponder network and 
the approximate position of Pegasus drop #157 positioned in a 


local coordinate system using Tl (12.5 kHz) as 0,0. 
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Figure 7. The C10 Local Coordinate System 


10 


cx DESCRIPTION OF THE PEGASUS SYSTEM 

The oceanographic instrument Pegasus is a free-falling, 
acoustically-tracked velocity profiler. It was developed by 
H.T. Rossby and D. Dorson at the University of Rhode Island a 
decade ago to fill a need for an instrument which would accu- 
rately record the fine scale vertical structure in the ocean 
and prove inexpensive and easy to handle at sea. It provides 
a means to measure absolute velocity components throughout the 
water column. (For the seminal paper on Pegasus, see Spain et 
aes, 1981.) Prior uses of Pegasus in studies of the Gulf 
Stream off the east coast of the U.S. have been documented in 
a paper and data report by Halkin et al. (1985). 

The position of the instrument as it descends and ascends 
is determined by travel time of signals emitted by Pegasus at 
10 kHz every 16 seconds which are heard and then answered at 
other frequencies by transponders previously set in place on 
the ocean floor, and by pressure readings. These time 
intervals, and temperature, conductivity, and pressure data 
from oceanographic sensors, are recorded and stored in a 
microprocessor-controlled memory in the instrument to be later 
down-loaded aboard ship when the Pegasus is recovered. 

The Pegasus model currently adapted for use by the Naval 
Postgraduate School is manufactured by Benthos. It is housed 
in a 17-inch glass sphere protected by a hard hat and, 
according to manufacturer claims (Benthos, 1989), provides 


operation to full ocean depth, integral flotation, no 
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corrosion, and large battery capacity allowing for 100 deploy- 
ments without opening. (A normal NPS school cruise uses about 
20 deployments.) It weighs approximately 45 kg in air and 
carries eight kg of expendable weight which sink the profiler 
in this present study at about 27 m/min. Six receiver 
channels for communication with transponders are available 
with frequencies located at 0.5 kHz intervals ranging from 
11. SekHZ "eo  14,,08knz. NPS supplied SEA BIRD model SBE-3 
temperature sensor and SEA BIRD model SBE-4 conductivity 
sensor are externally mounted with electrical interfaces, 
along with a  Paroscientific model 410KT pressure 
transducer with a companion Intelligent Transmitter circuit 


board. 


D. PEGASUS POSITIONING 
1. Computation of Pegasus Coordinates 
The position of Pegasus at each 16 second interval can 
be determined by solving a system of equations using the 


ex ap ihbUlrst 
Time = Distance/Sound Velocity 
Using one-way signal times from two transponders (see 
Figure 8), the depth of Pegasus as derived from pressure 
information, and an appropriate sound velocity, there will 


then be two equations in two unknowns, allowing for a solution 


eZ 
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Figure 8. Determining Pegasus Position 


of XP and YP at a given Pegasus position. The solution is 


derived from the equations: 


Time, = [(XP - X,)* + (YP - Y,)* + (2P - 2,)*))/7/v 
Time, = ({(XP - X,)* + (YP - Y,)* + (ZP - 2,)?))/*/V 
where: 
Time, = one-way travel time from transponders 1,2; 
V = effective sound velocity along the path between 
Pegasus and the transponder; 
XP,YP,Z2P = Pegasus position coordinates (ZP is assumed 
known from the simultaneous pressure 
observations) ; 


is) 


Nea = transponder coordinates (assumed known). 


Vipy! 


The traditional oceanographic method outlined above 
calculates Pegasus position coordinates. However, it lacks a 
means of obtaining precision estimates for the derived 
positions, and it is not able to use all the observed data 
Simultaneously. Furthermore, there iS no redundancy in the 
computational process and thus systematic errors or blunders 
can escape unnoticed. 

2. Estimating Horizontal Current Velocity 

Current velocity vectors u and v (describing velocity 
in the x and y directions, respectively) between two Pegasus 
positions, can be obtained by taking the difference between 


coordinates and dividing by the time interval: 


u = (XP, - XP,)/dt 


Vv = (YP, a= YP,) /at 


at = 16s 


This method, like the one mentioned in the previous 


section, does not provide estimates of precision and does not 


make use of all available data. 
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Lil. .SOURCESOF ERROR 


Davis et al. (1981, pp. 15-20) defines measurement error 
as the difference between a measured and a "true" value and 
describes these errors as being generally of three types: 

* blunders or mistakes. 
* systematic errors. 
* random errors. 

Blunders can be caused by carelessness, equipment failure, 
or false interpretation. They are usually large enough to be 
easily spotted when results are analyzed. 


Systematic errors follow a defined pattern or systen, 


which when discovered, can usually be expressed 
mathematically. Such factors as observer limitations, 
instrument imperfections Or inadequate calibration, 


meteorological conditions, or poor choice of mathematical 
model can produce the pattern which will remain consistent as 
long as the elements of the system remain the _ same. 
Systematic errors are not eliminated by repetition of 
measurements. They must be ferreted out and either removed 
from the observations or their effects added to the model. 
Thus, to begin such an analysis, serious consideration has to 
be given to locating and defining errors resulting from 


systematic and observational inaccuracies in the data. 


IS} 


Random errors are those that remain after mistakes and 
systematic errors have been accounted for. The values of 
these errors should be unbiased and will ideally distribute 
themselves about a mean of zero. It is these random errors 
that the least squares process attempts to minimize. In the 
analysis described in this thesis, a considerable amount of 


time was spent in determining what the estimates of the a pnori 


standard deviations of these random errors should be, both 
from the point of view of reasonable physical reality and 


meaningful results. The adjustment technique uses these a pnori 


estimates to weight the contribution of each observation to 
the final seiution. 

The following potential error sources were studied and, 
where applicable, appropriate standard deviations of the 
errors were then entered as weights into the least squares 
adjustment model: 

* Transponder coordinates. 
* Signal time. 
* Velocity of sound. 


* Pressure/Depth relationship. 


AG TRANSPONDER SURVEY 

Of prime importance is the precision of the transponder 
positions. Determining these coordinates is a difficult and 
time-consuming process. As noted before, in the C10 


transponder net, the location of the 4th transponder could not 
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be determined by traditional methods because the band width of 
the ship's receiver was not wide enough to pick up the 
ins KHz signal. Also, it was not possible at the time to 
determine more than two of the three baselines between the 
three "known" points, thereby precluding a mathematical 
closure of the triangle. 
1. Transponder Depth 

To establish a reasonable estimate for the standard 
deviation of transponder depth error, it was noted that any 
horizontal offset resulting from the ship's not being directly 
overhead always produces a positive depth error, i.e., the 
Slant range must be longer than the vertical. Figure 9 is 


Similar to Spain et al. (1981, p. 1557) and illustrates this 


process. 
The solution is as follows: 
(H + h)? = x? + H? 
nh ho = x 
a ar es 
where: 


H = true transponder depth; 


im 9 
Il 


depth error. 
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Figure 9. Transponder Depth Error 


If h << H, the last term can be disregarded, and 


= 
I? 
nlx 


A realistic estimate of the ability of a survey vessel 
to cruise directly over a transponder by using this method is 
about 100 meters (Schnebele, 1990a). Given the uneven canyon 
terrain of this survey which makes resolving the position of 
the shoalest depth even more difficult, it seemed reasonable 
to ascribe an error of 200 m to the offset. Thus, if the true 
depth lies at 2000 m and there is a horizontal offset of 


200 m, there will be a depth error of about 10 m. 
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Therefore, a o, = 10 meters was chosen for transpon- 
gers 1, 2, and 3% Since there was no observed depth for 
Transponder 4, a depth of 1900 meters was read from the latest 
NOAA chart using the position of transponder drop’ and, since 
it had not been surveyed in, a larger standard deviation of 


100 meters attached to it.? 


O(t1z,T22,T3z2) = LO mM 


Or, = 100 m 


2. Baselines 
In determining baselines, three sources of error 
predominate: 


* Minimum range reading taken when the ship is offset from 
the actual baseline. 


* Transponder depth error propagating into the computed 
baseline length. 


* Error in ship's course affecting baseline azimuth. 
In a simplified ideal situation where the baseline is 
crossed at a 90° angle in the center of the line and the 


depths of the two transponders are the same, the expected 


IThe transponder position was calculated from Loran 
Signals which are based on North American Datum 27. The NOAA 
chart uses NAD83. Positions may vary up to 100 meters between 
the two datums in this locality. 


“After the October 17, 1989 earthquake, a depth measure- 


ment of 1931 meters was made of the 4th transponder by NPS 
scientist Tarry Rago in a submersible. 
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error in baseline length from measured slant ranges can be 
calculated as follows (see Figures 10 and 11). 


a. Baseline Error Due to R 





Figure 10. Baseline Error 


R = distance off baseline when slant ranges were 
measured; 

B = true baseline between transponders (B = B, + By); 

jo’ = baseline error (b = b, + bz); 


Zo 
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This error increases the length of the baseline. 
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H = true transponder depth; 
ig! = depth error; 
S, = slant range from ship to transponders 1,2. 
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As shown in Figure 11, the hypotenuse representing 
the given slant range of the signal must be the same in both 
triangles ADE and ACF. FC must equal ED. Therefore, if the 
assumed depth is increased by error, then the computed 
baseline must consequently decrease. 

c. Total Baseline Error 


Thus the total baseline error bp + b, becomes: 
= il 2 
Drotal iy Dp + Di, ~ 52k 4Hh) 


Using typical figures of: B = 2000 m, H = 2000 n, 
h = 10m, and R= 100 m, it should be possible to obtain a 
value for by.4; = +30 m. 

If the transponder depths are not the same, and 
the baseline is not crossed right at the center but is crossed 


on a perpendicular heading, the baseline error will be: 





Protal pile 2°*B B 


In addition to baseline errors discussed above, an 
incorrect azimuth will add even more uncertainty. An azimuth 
error of 0.5° over a baseline of 2.0 km, for example, will 
produce an error in position of one end of the baseline with 


respect to the other of 18 meters. 
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Based on the above analysis, it was initially 
assumed that the relative horizontal coordinates of T1, T2, 
and T3 would be known to approximately +30 m. Hart (1967), 
for example, claims that with careful repetitive procedures 
involving all baselines in a transponder array, relative 
accuracies of the order of five meters can be achieved. 

As preliminary data analysis proceeded with the 
C10 network, it became obvious that there were major inconsis- 
tencies with the given baseline data. Preliminary adjustments 
(see Chapter V) would not converge when the given transponder 
coordinates were used. Further analysis suggested both length 
and azimuth problems with the given T2-T3 baseline. The fact 
that the T3-Tl baseline had never been determined led to a 
situation in which no positional redundancy existed in the T1, 
T2, and T3 network. As a result of these problems, initial 
positions of T1, T3, and T4 were determined from the Loran 
coordinates of their drop sites, while T2 was computed 
relative to T1 from baseline survey data. 

AS a consequence of this less than ideal 
procedure, it was considered that standard deviations of +100 
m should be allocated to the transponder coordinates of T2, 
T3, and T4, while T1 would be assumed fixed in order to 
provide a positional datum (albeit a rather arbitrary one) for 
the array. The 100 m allocated for positional standard 


deviations seemed reasonable in the light of both the known 
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positional accuracy of Loran and the added uncertainty from 
drift as the transponders settled to the bottom. 
In summary, therefore, the following standard 


deviations were assumed: 


B. TIMING 

The time observations required a consideration of how well 
the acoustic time intervals could be established between 
Pegasus and the transponders. The travel times were assumed 
to be independent of one another. The Pegasus manual states 
that the transponders have a 12.5 ms output pulse delay 
accurate to within +0.5 ms (Oceanographic Instrument Systems). 
For the initial transponder adjustment, o7;,. = 0.002 s was 
chosen as a reasonable weight (see Chapter V), and then later 
reduced to 0.0005 s for the full run with improved transponder 


coordinates. 


Ope OOS 


Ce OOO ams 
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Gee vELOCITTY OF SOUND 

Fofonoff (1963) states that velocity of sound is a 
function of the thermodynamical and chemical state of the 
water and is determined by using any complete set of variables 
of state such as temperature, pressure, and salinity. He goes 
on to observe that the precision of sound velocity must 
consider not only that of the oceanographic data, but also of 
the empirical formula used to convert these measurements to 
sound speed. 

Pegasus provides temperature and conductivity data which 
is converted to a sound velocity profile for the full range of 
the drop by algorithms published in the Unesco Technical 
Papers in Marine Science #44 (Fofonoff and Millard, 1983, pp. 
11-12, 49). 

In dealing with slant ranges, as in this case, acoustic 
refraction should be considered as well. Spain et al. (1981, 
pp. 1562-1564) presents an equation derived by Vaas (1964) 
which corrects for the effects of ray bending, taking into 
consideration the relative positions of the projector and the 
receiver even when they are at similar depths. This equation 
is stated to have an accuracy of better than 0.2 m/s for all 
angles and depths. 

Another excellent source for acoustic refraction informa- 
tion in a survey situation is the SASS Accuracy Study 
Simplified Ray Bending Correction (General Instrument Corp., 


1975) and its follow-up Ray Bending Correction for Depth 
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Sounders, An Informal Approach (General Instrument Corp., 
7G) 4 These two papers develop the equations for an 
effective sound velocity in detail and include graphs showing 
the errors relative to lateral angles. 

In this study it was found that the use of an average 
harmonic mean sound velocity produced acceptable results. To 
obtain this average, the sound velocity profile was used to 
calculate harmonic sound velocity profiles with respect to 
each transponder. These four harmonic mean profiles were then 
averaged (see Figure 12) and the resulting profile used to 
calculate the approximate X and Y coordinate of each Pegasus 
position. The harmonic mean profiles differ for each 
transponder because of the wide difference in deployment 
depths. Harmonic means at the depths found here differed by 
about 1.0 m/s from the average used, representing a bias of 
less than 0.7 m ina typical one km path length. 

Refraction was not regarded as significant for this study 
because of the relatively short path lengths and small sound 
velocity gradients encountered (Schnebele, 1990b). In deeper 
waters with longer paths, refraction may be significant. 

The harmonic mean used in this model, therefore, assumes 
no refraction. The small error that it introduces was felt to 
be inconsequential compared to the uncertainties in 
transponder coordinates and signal travel time measurements. 
Thus the harmonic means used in the adjustment were assumed to 


be exact quantities without significant random error. 
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Figure 12. Harmonic Sound Velocity 


D. PRESSURE/DEPTH RELATIONSHIP 

The algorithm used to convert pressure to depth is the 
Saunders and Fofonoff formula which uses the hydrostatic 
equation and the Knudsen-Ekman equation of state (Fofonoff and 
Peellard, (1983, p. 25). The formula includes variation of 


gravity with latitude and depth and assumes standard sea 
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water. This formula is said to deviate by only 0.08 meters at 
5000 decibars from estimates based on EOS80, a considerably 
smaller error than those in pressure measurements as shown 
below. 

The manufacturer of the Pegasus pressure transducer claims 
a typical accuracy under difficult environmental conditions of 
0.02% (Paroscientific, 1986). At our full scale of 2200 mn, 
this would compute to about 0.4 m. As the analysis 
progressed, further inconsistencies in the data began to lead 
to a suspected bias in pressure measurements, perhaps caused 
by a temperature hysteresis as Pegasus drops and rises through 
the deepest part of the water column. Therefore, a °p, = 2.00 


m was taken to be the standard deviation of the Pegasus depth 


measurements, i.e., 


30 


IV. LEAST SQUARES ADJUSTMENT 


The specific task of this study was essentially to analyze 
three basic problems: 

* The precision of the transponder coordinates, especially 

those of T4, from which the position of Pegasus is 


derived. 


* The precision of the Pegasus positions and consequently 
the current velocities derived from Pegasus navigation. 


* The depth at which pressure becomes critical in solving 
for velocities. 


To provide a means of answering these questions, a Least 
Squares Adjustment Program was developed by Dr. John Hannah.° 
This adjustment program, as well as a brief explanation, are 
provided in Appendix A. 

Least squares adjustment techniques are used to determine 
unique solutions for unknown parameters when there are 
redundant observations, i.e., more than necessary to specify 
the model (Davis et al., 1981, pp. 38-39). The least squares 
estimator is an unbiased minimum variance estimator which is 
unigue, mathematically easy to derive, and, when compared to 
other estimation techniques, leads to a smaller dispersion of 


Fandom errors. It also is distribution free in the sense 


“Adjunct Research Professor, CNOC Chair in Mapping, 
Charting, and Geodesy, U.S. Naval Postgraduate School, 
Monterey, CA., 1988-1990. 
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that a distribution is needed only for confidence interval 
testing (Hannah, 1989). 
The mathematical model for the Pegasus data adjustment has 


the general form: 


ct 
= 
| 


= Fy (X,) 


L.2 = F2(X,) 


in which the first set of observations comes from Pegasus 
itself in the form of return signal travel times from the 
transponders to the instrument. The second set comes from an 
a prio Knowledge of any of the system parameters such as the 
positions of the transponder coordinates, and each Pegasus 
depth as calculated from pressure. In general terms, any set 
of adjusted observations expressed as a function of a set of 


adjusted parameters can be written in a matrix expression: 


L, = F(x,) 


The following is adapted from a least squares adjustment 
procedure written by Hannah (1981) and is used here with the 


author's permission. 


...These functions may be linearized by taking a first order 
Taylor series expansion about some approximate values for 


ez 


the parameters, X. The first function then becomes 


] 
lye PV, =a | Cea) CXS) 
"a X_=X 
O 
OL 
A 
Ie + Vv; =_ A,X + |e 
Or 
A 
Vance Th, 
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oF) K 
Ar = sy 1» X = (Xa7Xo) 
¢ a xX io 
0 
and 


L, = Ly - 


Similarly, the second function may be linearized to give 


A 
V,.= A,X tg L» 


In the above, V, and V, are the vectors of residuals arising 
from the observations L,’ and L,* with the adjusted values of 
the parameters being given by the vector X,. Assuming that 
the two sets of observations are uncorrelated, then the 
least squares minimum variance estimate of X based on these 
two sets of observations is given by minimizing the function 


A A 
d = V,'P,V, + V"P,V, - 2K,'(V,-A,X-L,) - 2K,'(V2-A2X-L;) 


A 
Pienerespect to the wnknowns, Vv), V.,, K,, K,, and X. The 
weight matrices for each set of observations are given by P, 


33 


and P,, in which P, = 5-1 and Ps = x2 -1, i.e. the inverse of 
b b 


the variance-covariance matrices of the observation sets. 
Infinitely large variances are applied to non-weighted 
parameters resulting in zeros in the corresponding diagonal 
elements of P,. 


After minimizing and eliminating the unknowns K,, kK), 
Vi, and Vz, the least squares estimate for X is given by the 
solution of the normal equations 


A 
(A,'P,A, + A,'P.A,)X + (A,'P,|L, + A,7P.L,) = 0 


Since, however, the A, matrix arises from direct observa- 
tions on the parameters, the partial derivatives with 
ur 
respect to the parameters ax = I, the identity 
“a X =Xo 
matrix and thus the above equation reduces to the form 


A 
(A,'P|A, + P.)X + (A,'P,|L, + P.L,) = 0 
This has the solution 
A — 
X = -(A,;'P\AD + Po) FOCAL 


The variance covariance matrix of the parameter estimates is 
obtained by normal error propagation methods and is given by 


ae = (A,"P|A, + P,)7 


with the a posteriori variance of unit weight by 


gts T 
A Vy PSVy tVo PoV5 
Gy c Ne = 


in which n, equals the number of [Signal time interval] 
observations, n, equals the number of a priori parameter 
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observations, and u equals the number of parameters in the 
adjustment.... 


If the assumptions going into the adjustment are good, 


then the a posterior variance should be close to 1.0. This 


statistic is an indicator of the quality of the adjustment. 
The least squares adjustment technique is an excellent 
tool for determining solutions of unknown variables. However, 
certain limitations do apply. The process assumes that all 
systematic errors have been removed or accounted for, and that 
all remaining errors are randomly distributed. Systematic 
errors that still unaccountably exist in the observations will 
bias the adjustment such that it may attempt to distribute the 
errors to all the observations and shift the parameters. 
Other factors which may degrade the adjustment are: 


* A poor physical model such as one that ignores scale 
error. 


* The incorrect weighting of observations. 

* Small residuals which may be a result of poor network 
geometry or insufficient observations rather than good 
observations. 

For additional information on adjustment computations, see 
Uotila (1986) for derivations of appropriate expressions for 
least squares adjustments which use a variety of different 
systems or groups of observations with their associated 
constraints. 


Appendix A gives greater detail on the least squares 


adjustment process for this Pegasus data analysis. 
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V. COMPUTATIONAL PROCEDURE 


The data set used for most of the results below was the 
downcast data from Pegasus Drop #157 on August 3, 1989. 
Occasional comparison studies were also made with the #157 
upcast data, and that of #158 which occurred in the same 
vicinity later the same day. 

As stated in Chapter II, Pegasus data is recorded every 16 
seconds as the instrument descends and ascends during a 
deployment. For this study, the portion of the depth profile 
of Drop #157 from about 16 m to 2200 m was used, providing 306 
Pegasus records. This carried the instrument from surface 
waters down through the plane of the transponders. 

It needs to be stressed at the outset that the 
computational procedures adopted in this thesis should not be 
considered as optimum, but rather were developed as processing 
proceeded in order to overcome difficulties encountered with 
the specific data set associated with the C10 array. If the 
problems ultimately discovered with the data set had been 
known at the beginning, then some procedural points in the 
data processing would have been slightly revised. This aside, 
the purpose of the study was to demonstrate the capabilities 
of least squares estimation procedures to resolve both unknown 


transponder positions and Pegasus velocity components. As 
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will be seen in the following discussions, this was more than 
adequately achieved. 

Because of the poor initial positions of the transponders, 
it was decided to select a small data set of 74 Pegasus 
records from the total of 306 collected during drop #157 and 
use these to help provide improved Pranesondes coordinates. 
This data set was determined by dividing the full depth of the 
Pegasus drop into 36 intervals and taking a pair of 


consecutive records in each interval. Using the a prion 


standard deviations deeiea iM chapceremenl I vane initia | 
solution for the 74 Pegasus positions and the four 
transponders was obtained. The very small positional standard 
deviation associated with Tl essentially served to hold this 
transponder fixed in the resulting solution. 

As described in Chapter IV, (A,'P,A, + P,)+ is the variance- 
covariance matrix of the adjusted parameters. The upper 
tridiagonal portion of this matrix is stored in vector form 
and can be accessed to give information on the variances and 
covariances of the adjusted parameters as shown in Figure 13. 

For example, the portion of this figure which describes 
Pegasus Position 1 gives the variances of the x, y, and Zz 
coordinates, 0p,;*, Op7, and op,;*, respectively. The s other 


three elements in this upper tridiagonal, o and 


f e) a 
xl viene xl 21 


Op ‘lai provide the covariances between the x, y, and z 
vl Z 


coordinates of the first Pegasus position. Similarly, the 
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covariances between the coordinates of different positions are 
provided in their appropriate locations as_ shown. The 
variances and covariances of the transponder coordinates are 
located in the lower right hand corner, and among themselves 
comprise an upper tridiagonal 12 x 12 submatrix. 

The adjustment solution is also used to compute observa- 


tional residuals which in turn are used to compute an a posterion 
variance of unit weight. For the 74 record run this a posterior 


variance turned out to be 0.1196, a very small number compared 
to that which should ideally have been close to unity. This 


a posteriori variance when multiplied through the variance 


covariance matrix enables the matrix to be scaled in order 
that it provide supposedly unbiased estimates of the parameter 
variances and covariances. This was done with this data set 
and the resulting transponder coordinates with their newly 
estimated standard deviations (of the order of 10-20 m) taken 


and used as aprion information in the final solution in which 


all 306 records were used simultaneously. It was felt that 
this procedure would enable the transponder coordinates with 
their associated accuracy estimates to more closely represent 
the type of situation usually found in a good transponder 
network. 

In retrospect, however, it appears that it may have been 
more appropriate to have run the 74 point data set with oa;,,. = 


0.0005 s rather than the 0.002 s actually used. When, toward 


Bye, 


the very end of this study, a o,;,,, of 0.0005 s was allocated to 


the measured time delays in the 74 point data set, the a posteriori 


variance of unit weight became 0.96, and the resulting 
standard deviations on the transponder coordinates 
approximately 40m. The coordinate solutions for the 
transponders did not change by more than 5 m in position 
although the computed depth of T4 did increase by a further 
197 jie 

Although this second solution is to be preferred over the 
first, it was felt that the work and time involved in altering 
all the results and documentation already completed was not 
justified, given the minimal impact that it would have on the 
final results. In fact, it was clear that it would not change 
any of the conclusions resulting from this study. 

Ideally, when dealing with data from a number of Pegasus 
drops on a single transponder network, it would be best to use 
a small, but representative (in depth) data set from each 
drop, and merge these together in a single adjustment to 
provide an optimum set of transponder coordinates for that 
array which could, where appropriate, be held fixed in 


subsequent simultaneous processing of complete drops. 
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VO 


RESULTS OF STATION C10 ADJUSTMENT 
A. TRANSPONDER COORDINATES 


With improved transponder coordinates and their associated 


was then run. 


standard deviations determined through the 74 record procedure 
described in the last chapter, the entire 306 record data set 


The results from this 
discussed below. 


full analysis are 


2 below. 


The final transponder positions were moved (in total) from 
their original Loran estimates as shown in Figure 14 and Table 
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Figure 14. 


Transponder Movement 
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TABLE 2 


MOVEMENT OF TRANSPONDER POSITIONS 


Transponder X Y Z 

T1 0 0 4.5 
LZ =2.0.5 ~41.9 =6772 
T3 Be 2 a is ag @) — 7 
T4 eK OPES: 149.3 = Ono 


It must be stressed, however, that this data proved to be 
only internally consistent. Using the same transponder 
positions (as determined by the 74 Pegasus position data set) 
to adjust data from another drop #158 (which used only 13 
Pegasus positions), physically nearby and later the same day, 
produced somewhat different transponder coordinates although 
it converged within itself. The resulting Table 3 is shown 
below. 

It is suspected that the reason for this inconsistency 
lies both in the very low degrees of freedom in the adjustment 
of drop #158 (leading to a statistically weak solution) and in 


the overall weak aprior’ positions of the transponders. Ina 


well-surveyed four transponder array, this problem would not 
exist. 

It appears certain that the location of Pegasus within 
this weak transponder array has an effect on the final 


adjusted transponder positions. With the origin held fixed in 
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Drop #157 


Drop #158 


Difference 
Between 
#157-#158 


both cases, the array tended to skew in the direction of the 
drop (see Figure 15). 
appropriate for all drops might be obtained by making a series 
of Pegasus drops out along the edges of the array near the 


centers of the baselines, 


center. 


TABLE. 3 


teen oeONDERK COORDINATE SOLUTIONS 


trans— 
ponder 


T1 
T2 
T3 
T4 
T1 
T2 
T3 
T4 
T1 
T2 
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A set of coordinates which would be 
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Figure 15. Adjusted Array from #157 and #158 Drops 


It 1s well to note, however, how the transponder positions 
have been improved over those of our originally assumed 
coordinates, especially those of T4 whose position was largely 
unknown (see Table 4). 

The horizontal drms values were (when using ao; =0.0005s): 


‘ofa LBL 


O (held fixed). 
ee — Vee iN 
TS = 132 
x T4 = 12.4 m. 
These drms values suggest an improved horizontal precision 


of 13.5 meters or less. 
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TABLE 4 


STANDARD DEVIATION OF TRANSPONDER POSITIONS 


Transponder X Y Z 

Tl 0 0 10 
fai tial Tez 100 100 eG 
Standard 
Deviations ies 100 100 10 

T4 100 100 100 

a OO. O07 1.88 
Standard a2 8.32 7.48 75 
Deviations 
After 3 6.79 ios 27803 
Adjustment 

T4 11.67 Ate be 2./8 


Be PEGASUS POSITIONS 

The precision of the Pegasus positions is determined from 
the variance-covariance matrix as described in Chapters IV and 
V. Typical values at various depths are shown in Table 5. 

If geometric studies are desired, these variance- 
covariance values produce a tri-axial error ellipse for each 
Pegasus position. For ease of perception, the axes of these 
ellipses can be converted to an orthogonal system through the 
use of eigenvalues and eigenvectors. (For a discussion of 


this process, see Mikhail, 1976.) 
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TABLE 5 


PRECISION OF PEGASUS POSITIONS 


Depth (m) Oy Oy GO, 

16 Wee 632 O= Sal 
248 4.61 6.18 0.49 
503 4.31 6.00 0.56 
730 3.96 5.80 Oo 
1000 3%70 Se oe O58 
T250 3.44 SSS 0.67 
5 Oz Bias S16 a 0.86 
1750 3.24 5.26 WG Ae 
2005 B53 506 ee 
2183 3.80 4.76 I. oS 


Cc. CURRENT VELOCITIES 

Pegasus velocities were determined by using the adjusted 
values of the Pegasus positions which were determined by the 
methods described in Chapter V. 

It is well to note here the similarities and differences 
between the Pegasus positions as determined by this 
adjustment, and those determined by the method currently in 
use by the NPS Oceanography Department. 

The initial Pegasus positions used to start this adjust- 
ment had been derived by combining three observation equations 
into two, and then solving the two equations for the unknown 
x and y of the Pegasus position. (See procedure explained in 


Appendix A.) These positions were then refined by the 


46 


adjustment process using four transponders where Pegasus depth 
was constrained by pressure. 

The present method used by NPS as described in Chapter II, 
uses only two observation equations to solve for the positions 
(i.e., using only two transponders), without further 
adjustment and without the pressure/depth constraint. 
Velocities derived from these two sets of positions are 
compared in Figures 16 and 17. (Both sets of velocities have 
been filtered by a seven point running average.) Figure 16 
plots the U components of adjustment derived velocities 
against depths, those derived by using the present NPS method, 
and finally the differences between the two approaches. 
Similarly, Figure 17 plots comparable information for the V 
components. 

The profiles of the two methods exhibit very similar 
configurations through much of the range, with notable 
exceptions occurring at the surface, at mid-range (see 
discussion of residuals in the section on depth/pressure 
analysis which follows), and at the bottom of the cast. 

The differences in surface velocities between the two 
methods may be a result of using different sound velocity 
profiles. The large divergences at depth where the adjusted 
speeds peak at about 22 cm/s and those of NPS at 111 cm/s, are 
almost certainly due to a loss of precision in the two- 


transponder solution, as well as inconsistencies in the 
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Figure 16. U Component 
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Figure 17. V Component 
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Original pressure measurements. Table 6 provides a sample of 
the velocities derived by the two approaches, and shows the 
differences between them. These differences begin to increase 
markedly at about 2100 m, indicating the substantial role the 
adjusted transponder network and the inclusion of a pressure 


observation play in the solution. 


TABLE 6 


VELOCITY COMPARISONS 


Adjustment Velocity Differ- Differ- 
derived NPS Differ- ence ence 
Depth velocity* Velocity rence in U in V 
(Approx. ) (cm/s) (cm/s) (cm/s) (cm/s) (cm/s) 
m 1 2 Ze )een (elk 
40 14.23 E3324 = 10) SNS, +630 +5.49 
P12 oe. oo 25.41 +652 he. OO +6222 
660 210 OS i. 45 =0:..35 +QO72 2 
1280 10.01 10.74 +0.74 pale 0.58 =O 77 
1380 1.08 9.45 TG. =7 706 =] 2 
1820 10.34 £24166 WZ «53 tise oo +1.30 
2090 £7.32 aoe 2 +2416 1 OnG.1 +2 20 
2112 19.58 35.84 Gree) sad Co betsy Jl +11.44 
2145 23707 68.05 +44.98 +3328 9 +30.10 
2175 6236 110.75 O45 oO t+ 71.56 +6252 
2190 Mess Sis) 49.76 +33.40 tee 96) +22 562 


*Adjustment using o, = 0.002, 
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ee 
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Seven point filter used for both sets of velocities. 


D. PEGASUS VELOCITY ERRORS 

Estimated velocity errors were obtained through the 
propagation of positional errors found in the variance- 
covariance matrix of the final run, yielding a,, o,, and a,, the 
errors in the x, y, and z directions, respectively. This 
propagation of errors proceeds as follows. 

The desired velocity components are given by the 


expression: 


7 U (Gea) 7 oe 
V= |Vj] = (Yo-yi) At 
W (Z.-Z,) /At 


where (X,,Y1,2,;) and (X2,Y2,2Z2) are the positions of Pegasus at 


respective 16 second intervals (i.e., At = 16s). 


The variance-covariance of V is given by 


Ou Guy O uw 
ee 2 2 i 
Lz = Ov, Oy se = G2Z,pG 
Gem - (Cane 


where 


3/ OX, 0/ oy, 0/325 3/ 9X. 3/ 0Y5 3/92. 


= 0 0 | 


if 
a= | o - 
At zi 
0 0 -1 O 


oye 


and 


In this last expression XP) and XP, are the variance- 


covariance matrices for PegasuS positions 1 and 2, 


respectively, while XP, P., and 2 PoP, are their cross 
covariances. These are obtained from the full adjustment 
variance-covariance matrix described in Chapter V. (Appendix 


B contains an algorithm for computing the variance-covariance 
Matrix of the current velocities from data obtained from the 
subroutine MAT in the Least Squares Adjustment program.) 

These values provide information on the precision of 
Pegasus velocity at each position, and change according to the 
geometry between Pegasus and the transponders. Table 7 
provides velocities and their standard deviations at various 
depths as computed by the adjustment. The horizontal velocity 
errors are greatest near the surface, and improve as Pegasus 
descends toward the transponders where, in the horizontal 
sense, the geometry becomes tighter. For the vertical 
velocity, the reverse is true. 

These velocity errors were disappointing. At the surface, 
a velocity in the ‘*X direction of -17.8 cm/s had a standard 


error of +10.7 cm/s. In the Y direction, a velocity of 


bz 


TABLE 7 


VELOCITIES AND STANDARD DEVIATIONS 


Av. 
Depth 
(m) U (m/s) V (m/s) W (m/S) 9, (m/s) Go, (M/S) 9, (M/S) 
20 -0.1784 0.0169 0.4347 OREO 73 0.0822 0.0403 
Zo20 —0. 0886 0.0789 0.4566 0.0988 0.0759 0.0408 
740 Or. 2 316 0.0204 0.4511 ©2707 95 0.0621 0.045 
M733 -0.0729 0.0011 0.4141 0.0643 Of05 1 i 0.0573 
i 39 Orel C—O. O02 26 0.4734 Of 0575 0.0439 O21 107 
2180 0.1989 0.0681 Gris 379 O72 0569 0.0454 O23 359 


1.7 cm/s had a standard error of +8.2 cm/s. At the bottom 


results were somewhat improved with u = 19.9 cm/s and oa, = 5.7 


cm/s, v = 6.8 cm/S and o, = 4.5 cm/s, but not substantially so. 


These velocities and error estimates were obtained by 


looking at only individual pairs of positions. The NPS 


Oceanography Department computes velocities by using a seven- 


point running average for each Pegasus recorded depth. This 


same filtering technique used with the adjusted data would 


improve precision by 1/7 or 0.38, if positions were 


uncorrelated. However, these positions are correlated, so the 


expected improvement is somewhat smaller. 


There are two ways the precision of the unfiltered 


velocities could be improved: 


so that as At becomes larger 
the 


* Use longer time intervals, 
in the error propagation equation discussed above, 


variance-covariance matrix Z V diminishes. 


D3 


* Average the 16-second velocities (in the same manner as 
NPS) by using a seven-point running average which would 
improve the precision Lf the velocities were 
uncorrelated. 

However, both of these techniques, while gaining precision, 
sacrifice vertical resolution. 

There are two important aspects to these formal velocity 
errors. In the first instance the signal time intervals are 
only measured to 0.0001 s. At a conventional sound velocity 
of 1500 m/s this is equivalent to 15 cm in the derived range. 
This in turn propagates into a velocity error of approximately 
1 cm/sec between consecutive Pegasus positions if the 
transponder positions are assumed to be without error. 

In the second instance, the formal errors in the 
transponder positions propagate directly into errors in the 
Pegasus positions and thence into the velocity components. 
These latter errors have by far the most significant influence 
on the derived standard deviations of the velocity components 


for Pegasus. This in turn provides an added emphasis on the 


need for a strong transponder survey. 


Ee DEPTH/PRESSURE ANALYSIS 

To determine at what depth the pressure/depth measurement 
becomes critical in solving for current velocities, a separate 
computer run was made with the standard deviation of the 
Pegasus Z coordinate changed from 2.00 to 50.00 nm. The 


results are graphically illustrated in Figure 18. 
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Here, the U and V velocity components derived from the 
adjustment using op = 2m are plotted on the left; speed 
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differences between the U and V components of the two sets of 
adjustments using Op, = 2, and op, = 50, respectively, are 
lllustrated on the right. The differences are minimal down 
through the water column until about 1700 m as Pegasus 
approaches the plane of the transponder. DU reaches a maximum 
positive difference of 1.16 cm/s at 2008 m, and subsequent 


Maximum negative difference of -3.53 cm/s at 2135 m. 


ahs. 


Respective values for DV are 0.49 cm/s at 2001 m and 
—-2.26 EM/7Seat .245 5 on. 

An analysis comparing the standard deviations of Pegasus 
depth positions shows a large maximum of 36 m at a depth of 
2000 m as derived from the free-floating P, adjustment, 
compared with a maximum of 2 m at a depth of around 1940 m for 


tightly-held P,. The latter reflects the fact that the a priori 


precision estimate for pressure observations was itself 2 mn, 
and that travel time measurement geometry gives no additional 
information on the adjusted z value. 

Table 8 shows observational residuals after both of the 
adjustments and compares the two sets of residuals at similar 
depths. (These figures compare runs made with o,;,,, = 0.002 s, 
rather than the final 0.0005 s.) 

At the surface both sets of residuals are small. Where 
the pressure/depth relationship is tightly constrained, 
residuals were highest at the bottom of the cast where Pegasus 
moved through the plane of the transponders, i.e., between 
about 1870 m to 2240 m. Here the resolution of depth becomes 
less certain due to the poorer geometry of the Pegasus- 
Transponder array. These large residuals at depth indicate 
there may be a bias in the pressure measurements, thus 
suggesting the presence of a systematic error in the pressure 
sensor (perhaps due to a hysteresis in temperature readings) 


that is unaccounted for in the adjustment model. 
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Z Op. TT 2 Ts T4 
Surface 2 -0.0001 -0.0002 = (0. (OKO 0.0001 
0.0000 0.0000 0.0000 0.0000 
0.0000 0 (QOL -0.0001 -0.0001 
50 0.0000 -0.0002 0 0 
0.0000 0.0000 0 O. 
0.0000 Oe Oman 0 0. 
1400 m 2 -0.0004 -0.0012 Cie 0 
(approx. ) -0.0005 -0.0013 On 0 
-0.0005 -0.0012 On 0 
50 -0.0004 -0.0011 e. 0. 
-0.0004 -0.0012 0. 0 
-0.0004 -0.0013 0. 0. 
Bottom 2 0.0014 0.0011 0 0 
0.0016 0.0014 0 0 
0.0019 0.0016 On 0 
50 0.0001 0 0 
0.0002 0.0002 0.0000 -0.0001 
0.0003 -0 


Residuals at the bottom of the cast from the adjustment 
where depth was allowed a very loose constraint are uniformly 
low. The differences between these two sets of residuals at 
the bottom, as well as differences in velocities and an 
increase in — at this depth as discussed above, all tend to 


corroborate the suggestion of a bias in the pressure 


measurements. 


Sy) 


Residuals also increased in both cases for a few positions 
around the depth of 1400 meters. This appeared to occur where 
the trajectory of the instrument reached its most westerly 
position. Reasons for this are unclear, although a possible 
explanation may be related to the fact that this depth and 
position approximate the location of the shelf edge of the 
canyon. 

It is also possible that one or more of the travel time 
observations in this part of the drop had large errors. This 
interpretation is further suggested by the anomalous current 
shear computed at this depth from the original data (see 
Figures 16 and 17). 

It should additionally be noted that these large residuals 
at depth and at mid-range seem to display a systematic 
structure wherein their sign is consistently the same. If the 
errors were truly random, the positive and negative signs of 
the residuals would be distributed more or less evenly. This 
Situation is another indication of some systematic error which 


has not been accounted for in the adjustment model. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


It is unfortunate that the data set and the transponder 
array used in this study lacked the quality desirable for de- 
tailed analysis. Many of these difficulties were, unfortu- 
nately, only discovered as data processing proceeded. It is 
encouraging to note, however, that the least Squares proce- 
Gures described here enabled the identification of problems 
with the data which would otherwise have passed undetected. 

From a theoretical standpoint, the least squares process 
must provide better results than the standard NPS positioning 
techniques because of the fact that it uses all the observa- 
tional data. It has the added advantage of providing vital 
statistical information on the quality of the derived results. 
Unfortunately, in the case of the C10 network, the uncertainty 
regarding the original positions of the transponders made it 
Gifficult to perform a comprehensive comparison between the 
least squares technique and the standard NPS methodology. 

It is observed, however, that the least squares method 
obtained a solution for transponder positions that converged 
to less than one meter, with standard deviations that were 
much smaller than those with which the adjustment started. 
Precision of the transponder coordinates showed horizontal 
arms values of less than 15 meters, with standard deviations 


of transponder depth less than three meters. While these 
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figures may be optimistic (see Chapter IV), they are 
indicative of the improvements which can be achieved using 
these methods. 

Regarding the precision of the horizontal Pegasus 
positions, the standard deviation of the X values ranged from 
4.9 down to 3.2 m, reaching a minimum around a depth of 1800 m 
just at the upper limit of the transponder plane, and then 
increasing slightly to the bottom of the cast. Pegasus Y 
values showed standard deviations of 6.3 to 4.8 m, with 
passage through the transponder plane not as noticeable. 

Standard deviations of the Pegasus Z values, when they 
were tied tightly to pressure, ranged from 1.4 m near the 
surface to a maximum of 2.0 m at a depth of 1943 m. When the 
depth was only loosely constrained, standard deviations varied 
from 2.3 m near the surface, to 8.6 around depths of 1800 m 
just prior to entering the plane of the transponders, and 
reached a maximum of 36.5 at a depth of 2000 m. From there 
they steadily diminished back to 10.4 at the final depth of 
2165] 

Comparison of the results of the NPS method and the least 
squares method shows a considerable difference in velocities 
at depths where Pegasus passes through the plane of the 
transponders. This difference at depth occurs also in 
velocity comparisons between adjustments made holding depth 
constrained with a standard deviation of 2.0 m, and allowing 


it to float more freely with a standard deviation of 50.0 m. 
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This leads to a suspicion that there may be a systematic error 
in the pressure reading that has not been accounted for in the 
model, perhaps due to a hysteresis in temperature readings. 
Accuracy of pressure observations becomes critical at depths 
below 1700 m. 

It was disappointing to note that the standard deviations 
on the computed Pegasus velocities (when using At = 16 sec) 
were often of a similar magnitude as the velocities themselves 
(i.e., 5-10 cm/s). However, with resolution of signal travel 
time possible only to the nearest 0.0001 s, resolution of 
range becomes +15 cm. This in turn propagates into velocity 
errors in the order of 1-2 cm/sec. Therefore, it is unrealis- 
tic to look for much greater precision than that. 

It is recommended that: 

* NPS refine the existing least squares techniques and the 
adjustment software to facilitate its use on a regular 
basis, for production operations. 

* Four transponders be used in each network. This will 
both strengthen the solutions for the Pegasus velocities 
and provide a reasonable measure of redundancy for the 
adjustment procedure. 

* More acceptable positioning of transponders be under- 
taken, within the time constraints involved, including 
observation of all baselines. Close attention should be 
paid to both length and azimuth. 

* That a small, but well-distributed data set of Pegasus 
records be used from each drop to assist in providing a 
unique set of transponder positions for each network and 
that these positions be held fixed in the subsequent 


adjustment of the full set of records from each drop. 


* Where practical, pressure information always’ be 
collected. It both adds to the overall redundancy in the 
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* 


network and enables additional parameters to be 
introduced if necessary. 


A very well-controlled set of experimental data be 
collected on a well-positioned five transponder network. 
This will enable clarification of the following issues: 


- The lack of consistency between the adjusted transpon- 
der positions when using different data sets for the 
Same area. 


- The question of pressure/depth relationships, 
especially at depth. 


- Possible causes for high residuals at a given depth 
(1400 meters in the case of station C10). 


- The introduction of additional parameters to describe 
possible mis-calibration of the pressure head or 
hysteresis of temperature readings. 


- The investigation of the actual motion of Pegasus as 
1t descends and ascends. 


- The geometric impact of drops made on the edges of the 


transponder array aS opposed to those made at the 
center of the array. 
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APPENDIX A 


LEAST SQUARES ADJUSTMENT PROGRAM 


The least squares adjustment equation for a combined 


system of two sets of observations is (Uotila, 1986, p. 97): 
ss T T 7a i 
ee Cee ety tA AD) (AG Py + AZ P21) 


where, in the terminology of the program: 


A 

x = column vector of corrections to the parameters 
(NP, 1) j 

1 = Observations (NOBS,1); 

ix = Jacobian matrix, or the partial derivatives of each 


observation with respect to each parameter 
(NOBS,NP) ; 


Ae = Transpose of A (RP’,NOBS) ; 
P = Weight matrix (NOBS,NOBS) ; 
N = Maximum # of Pegasus positions to be determined; 
NT = Maximum # of transponders in the network; 
NOBS = Maximum # of observations 
NP = Maximum # of parameters allowed = (N x 3) + 
(NT x 3). 
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A. FIRST SE? OF O8SER ALIONs 


1. DL, Matrix 


The first set of observations uses one-way travel 
times from the transponders to Pegasus as the observations. 


Parameters are the x, y, and z coordinates of each Pegasus 


position. 
L, = F,(x) 
Time = Distance/Velocity 
Position 1: 
Time, = [ (XP,-X,) ° 7 Obits) a (ZP,-24) 7] /Vv 
Time, = [(XP,-X,)* + (YP,-Y2)* + (ZP,-22)2]?/V 
Time; = (XPj=x,) ob (YP i pee (ZP,-25) 2] ?/V 


iL 
Time, = {(XP,-X,)* + (YP,/3¥,)2 + (ZP,-Z,)7] 7¥ 


Position N: 


Time, = time intervals from transponders 1,4; 


V 


sound velocity; 
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XP,YP,ZP = Pegasus position coordinates; 


Mee, Dy = transponder coordinates; 
# Observations = N x 4. 
Each Pegasus position therefore has four time-interval 
observations, one for each transponder. Thirteen Pegasus 


positions will provide 52 time observations and (13 ~x 3) 
parameters, and therefore 13 degrees of freedom for this first 
set of observation equations. For each additional Pegasus 
position added to the data set, an additional degree of 


freedom is gained. If apron constraints are introduced for 


the transponder coordinates and the transponder positions 
solved for in the adjustment (as has been done here), then the 
number of observations and the number of unknowns rises by 12. 


The adjustment procedure requires an a prion knowledge 


of the transponder positions, the depth of each Pegasus 
position (derived from the pressure information) and the sound 
velocity for each Pegasus position. Estimates of the x and y 
coordinates are obtained by taking the three observation 
equations for each Pegasus position that contain information 
from transponders Tl, T2, T3 (presumably the best known). 
Subtracting the first equation from the second, and the second 
from the third, will yield two equations in two unknowns, and 
therefore a resulting first guess solution for XP and YP at 


that position. 
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The L, matrix in the adjustment equation is a matrix 
of differences between observed time intervals and what the 
observation equations would yield using our first guess 


parameters. 


Ly = L, calculated ~ L, observed 


2. eee eters on 


The A, matrix contains the derivatives of each 


observation with regard to each parameter: 


A, Matrix = cee 
OX 


i 


T= ( (XPj=X))° + (YPJ=Y))7 4 (ZP,-2,)2]?2/V — Vi 


i =A aaa 
ea 1 eee 
Then 
2B 91 BLO 
IXP. = 11 
del] DV 
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The derivatives of the transponder positions will be 


the negatives of the corresponding Pegasus position 
eerivatives. 
oF. (“pak ..) 
pol DV 
a3 
This Will ecreate a Matrix.1n. the form: 
e/a 0/0 
o/ >a Ye er, a/ dT) 3/ aT, a/ dT, a/ oT , 

ty Kot, OX x,x,xX 
Peg. Lo x, xs ee x 
ros ta~|/x,x,x x,xX.x 

ty Kee Ke ox 

ty Gre. Gen x, xy xX 
Peg. to eg Xs 5G 
Pos ts x,xX,xX xX ,xX,xX 

Uy x ,xX,xX Keen oN 

Sy Me Mik OCR e 
Peg. to x. ex ore GPP 4 
ros N ts x,xX,xX Ke Xa 

Ly x,xX,X Ke kx 
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where: 


t, = Time,; 

DP; = XP) Yee 
TT, = Xi, Yu, 27 

7 = 2 ee a: 

5 = to Ne 


If submatrices are described by 


Ke Go xX,xX,xX;,0;0,0) 0, O40 0 O70 
Bly Ke oe Sie = 0,0;0,x%,%,;%,0,0;,0,0,4070 
are en 4 0,0,0,0,0,0,x%,x%, 4,070 70 
» Pe ee, 0,0,0,0,0,0,0,0,0,% = 
Then the A, matrix becomes 
ai 0 0 e+, You > 1 oes 0 Si 
0 A> 0 0 So 
0 0 a3 0 S3 
e S e e e e & e e e e e e e e An Sy 


For each Pegasus position there will be a 4 x 3 block 
for the Pegasus derivatives and a 4x 12 block for the 
transponder derivatives. All other positions in the matrix 


will be zero. Therefore, to drastically reduce computer 
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memory and time requirements, the program was developed to run 
in a sparse matrix form whereby only the non-zero elements are 
used. To do this, the whole A matrix is divided into parts: 
* A--which contains the values for the Pegasus components. 
* S--which contains those for the transponders. 
3. P, Matrix 
P,, the weight matrix, is set up as a diagonal matrix 
and stored in column vector form, assuming no covariances 
between observations. If it is desirable to add covariances, 
then an error propagation subroutine can be inserted to fill 
out the weight matrix and the program altered accordingly. In 
our case the chosen standard deviation for the time interval 
was 0.002 s, later changed to 0.0005 s. It was entered into 
the program in the error propagation section where the code 
reads: 
Do 35 I =1,NT 
P(1,1) = 1.0D0/0.0005D0**2 
35 Continue 
B- SECOND SET OF OBSERVATIONS 
1. L, Matrix 
The second set of observations corresponds to various 
parameters which have been observed, in this case the position 
coordinates of the four transponders, and the position of the 
Z coordinate of each Pegasus position which allows the depth/ 


pressure relationship to be constrained. 
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Ce eee 
2S Fas 
i; =a 
x X,Y,Z of each transponder 
Lig = ea 
ZP, = ZP, Depth of Pegasus as computed 
1 = 1,N from pressure. 


As described above, we already have first guess 
estimates of all these observations. The L, matrix will 
initially be zero since for the first solution, the observa- 
tions equal the parameters, i1.e., xX, = X;. However, with 
further iterations, this will not be the case. 

2. A, Matrix 
The A, matrix contains the derivatives of each 


observation with regard to each parameter: 


For example: 


Jall else 


78, 


This produces the following diagonal matrix of 1's and 


O's 
O 
e) 
i 
0 
O 
1 
1 
1 
1 
ch 
1 
a 
ik 
alk 
ah 
1 
ie 
i 
Beeb) Mati x 
P, 1S a Similar diagonal matrix stored in column vector 
erm . It weights all Pegasus z positions to allow for 


tightening or loosening the depth/pressure relationship. This 
chosen standard deviation, in our case a value of 2.00, enters 
the program through an interactive request at the beginning of 


the program run. The standard deviation of each transponder 


var 


position coordinate is read into the program from the original 


data set. The P; matrix will @ppear as: 
0 
0 

] 
oe 
Zi 
‘a 

0 


if, 2 
Toe 


val 


Ye 2 


24 


C. MATRIX OPERATIONS 
Now all the matrices are in place and matrix operations 


can begin. 


Oe T T -} T T 
x (A, P\A, + A, P2A2) (A, P,L, + Az P22) 


(ee 


Since the A, matrix is a diagonal of O's and 1's, the 


second term in the first parenthesis merely adds weight to the 


appropriate diagonal position in the first term. 


D. REQUIREMENTS FOR RUNNING THE PROGRAM 


1. Initial Steps 


Before the program is run, the following steps must be 


taken: 


* 


* 


Variables must be dimensioned adequately, following the 
definitions clearly stated in the preface to the program. 


The value of o,,,. must be defined (in our case 0.0005), 
and/or the appropriate error propagation subroutine 
added. 


A data set file must be available for access and set up 
in the following fashion: 


- Transponder coordinates, each with its own standard 
deviation. 


- For each Pegasus position: pressure, depth, four 


transponder one-way time intervals, and sound 
velocity. (Sample data set is shown below.) 
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Transponder | Tl 
Coordinates |T2 
14 
14 
Transponder |T 1 


Standard 2 
Deviations {73 
T4 


oo Oo Oo 
- oO O89 oO ODO OO 


Peoasus ) 
Positions 2 
3 


fy 


. 


3 
iG 
.0 

iS 
ae 
.0 
8 
A) 
Fé 
D 


HW HR Tew WM Hh Mwon HN NN OH 
Sooo. OF CFC. Cee ee eS ae 


~ ND ORD 





T i Sound 


Press- | Depth Ua | di 
Velocity 


Ce 





Time Intervals 


* FILEDEFS for input and output files must be in place, 
either typed in before the program starts or available 
through a program exec. The program will run in either 
WF77 or FORTVS2 fortran. 

- FILEDEF 01 DISK fn ft fm (1i.e., Test Data A). 
- FILEDEF 02 DISK(recfm vb lrecl 132 blksize 134). 


- (For large data sets, use FILEDEF 02 DISK(recfm fb 
lrecl 132 biksize-13200)%- 


* Be sure enough memory is available. Two M sufficed for 
running at least up through 74 Pegasus positions, but 
4096 K were required for the 306 position run. 


The program starts off by asking for operator input 


(sample answers in parentheses), requesting: 
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* Project or run description. 
* Number of transponders (4). 


* Global standard deviation to be given to the Pegasus 2 
coordinates (2.00). 


* Convergence limit for adjustment and number of iterations 
(entered separately) (1.000), (4). 


From then on the program works on its own, finally 
producing whatever output has been requested in the output 
file. 

2. Brief Program Outline 
* Matrices are zeroed out at appropriate times. 
* Pegasus data set is read in. 


* Subroutine APPRO is called to calculate initial approxi- 
mate coordinates of Pegasus. 


* The diagonal weight matrix P, is set up and stored in 
matrix BIGP. 


* Program constants are calculated. 
* Weight matrix P, 1S generated. 
* The L, matrix is formed. 


* The matrices A and S are formed by calling subroutine 
Serer. 


* The normal matrix equations are formed block by block by 
calling subroutine NORM. Outputs include the normal 
matrix (ATPA), in column vector form called anorm, and 
xhat the ATPL vector called xhat (which subsequently 
becomes the solution vector). 


* Subroutine P2L2 is called to compute the P,L, matrix and 
add it to the column vector. It also increments the 
diagonal elements of the normal matrix to include the 
influence of P2. 


* The normal equations are solved by calling a user 
Supplied routine that will accomplish the required matrix 
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operations for producing the adjusted values with which 
to correct the parameters. The IMSL routine LINV3P was 
used in this study. Utilizing a symmetric storage mode, 
this algorithm replaces the ATPA matrix by its inverse 
which is the variance-covariance matrix of the Pegasus 
and transponder positions. This (ATPA)™ matrix can be 
written to file for later access in order to produce the 
error matrix required for current velocity analysis (see 
Chapter V) as used in the program provided in Appendix B. 

* The parameter corrections are tested to see if iteration 
1s required. If so, the program then updates the 
parameters and goes through the process again for as many 
iterations as are called for, or until the convergence 
limit has been reached. 


* Residuals and the apostenon variance of unit weight are 
computed by calling subroutine RESID. (Residuals here 
are, of course, given in units of seconds since they 
relate to time interval observations.) 


3. Program Outputs 
Outputs include: 


* Original data set and interactive information given at 
the beginning of the program. 


* Adjusted parameters following each interaction. 

* Residuals and a posteriori variance of unit weight. 
* Any other information requested to be written to a file. 
For this study a piece of code, Subroutine MAT, was 
attached at the end of the program to pick out from the anorm 
vector only the data required for current velocity analysis. 


This subroutine is called after the write statement for the a 
posteriori variance of unit weight which uses format statement 


7 Ly Se The anorm vector is the upper tridiagonal of the 
(ATPA) =? variance-covariance matrix stored in column vector 


form. (See the diagram in Chapter V.) 
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The position of the variance of each parameter can be 


obtained by: 


N x (N +1) 
2 
where N is the place number of the parameter. For example, 


suppose the variance for the Pegasus y coordinate in the 
second record was requested. YP, is the 5th parameter--xXP,, 
Peer, XP>, YP>, .--, (5*6)/2 = 15. The variance for YP,, 


therefore, will occupy the 15th place in the anorm vector. 


XP, YP, Ze xP, ves | ZP, 


di 2 4 7 11k 16 ee : 
Peg. Pos. 1 3 5 8 2 17 See 
6 9 13 18 ck = 
10 4 INS. 
Peg. Pos. 2 as) 20 5 5 € 
24. 5 Ane 


The subroutine Mat picks up only those values in the 
anorm vector that fill the 3 x 3 variance-covariance matrix of 
each Pegasus parameter, and the 3 x 3 covariance portion that 
relates each parameter to the next one in line. A sample of 


the output follows: 


Te 


24.2487314359615 3942 18.6328080418266140 0.934616249631990276E-01 


Var-covalr 18.63280804618266140 $9.8923498357777184 -0.591505777505010405 
mos 0.9346162469651990276E-01 -0.591505777505010405 0.260716185532883254 
ec 22.7007923927885546 18.5979145727731243 0.6312602313846052462E-0) 
abies 18.5612829360100449 3$9.0364503165567890. ~0.4466420607799607621 
Pos 1-Pos 2 0.85188 766419506B80S2E-01 -0.45395290611846188850 0.520710922203657129E-01 
Var-covar 24.099975088463 7862 18.609687082718410958 0.7239174621678468885E-01 
Pos ? 18.60987827184610938 39. 9102237605 952366 —-0.595096665108637957 
0.723917462167846885E-01 -0.5950966651086537957 0.260009015573909333 
c 22.567175 71932801941 18.570915890451900) 0.485779418760951404E-0} 
aes 16.550950735 9568424 3$9.04694978092177121 ~-0.445954216385781299 
Pos 2-Pos 5 0.6S5S6726100041890147E-02 =-0.4434994221539607538 0.510857467937755941E-01 
Var-covar 235.96487646886889567 18 .594349906284654} 0.565971236676837824E-01 
Pos 3 18 .5943499042846341 39.909825695841446) ~0.595998618935371994 
O.S6S97I25S66768S7B26E-0} -0.593996618955571994 0.259193690776075099 


The first set of three lines with three values each is the 
Variance-covariance matrix for Pegasus position 1, which is 
symmetric. The second three lines of three values is the 
covariance matrix for position 1 and position 2 
(nonsymmetric). The third set of three lines is the variance- 
covariance matrix for Pegasus position 2 (symmetric), and so 
on. 

The final 30 lines in the output provide the variances 
of the transponder coordinates and all the covariances between 


them. 
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Girne) Cooma C2 ©) 62°C) ©) ©) GC) Gat 0) G0) (2 0) C2 Co ta 02 ae kp 1 > I Ge Se | 


‘?) 


Gy @-Cic) 


PROGRAM PEGASUS 


ToarcorvYoicookccesdk tcc kak kok kote ak teat ak teat ak kak ak ak te ak te teat seat ak te ak ae ae ae ake aie re ae He aie ae ak ake ae ae ae ae 


at 
re 
* 
* 


DADADDDODADDADADDADODDANDADDDDAM® 


COMMENTS GIVING A DESCRIPTION OF THE PROGRAM ¥ 


Mokke kkk teak oka tk tolokack fever ake teat kak teak ke te fete te kk ake ate ate ae ak aie aie ak ae aie sk ae He He ke Neve ake ake Here Ye 


INPinieT REAE*8 CA=H,0-Z) 


Pe sooo eee 


SET MATRIX VARIABLES TO THEIR MAXIMUM SIZES. THE VARIABLES 
USED ARE AS FOLLOWS: 


N = MAX. # OF PEGASUS POSITIONS TO BE DETERMINED (20 ) 
NT = MAX. # OF TRANSPONDERS IN THE NETWORK (4 
NP = MAX. ## OF PARAMETERS ALLOWED = (N*3) + (NT*3) 
NOBS= MAX. # OF OBSERVATIONS = N*NT 

NV = MAX. # OF NON-ZERO VALUES IN THE A MATRIX = 6*NT*N 
NN = MAX. # OF ELEMENTS IN THE UPPER TRIANGULAR PORTION 


OF THE NORMAL MATRIX = NP*(NP + 1)/2 


THE FOLLOWING MATRICES MUST BE DIMENSIONED TO THEIR MAX. 
SIZE PRIOR TO THE PROGRAM BEING USED. 


XO(NP), P(NT,NT), BIGP(NT,NOBS) ,LB(NOBS,1), TRANS(NT), 
X(NT) ,YCNT) ,Z(NT) , ACNT, 3) , SCNT, NT*"3) , VALUE(NV) , ATP(NT, 3), 
ATPA(3,3), ATPS(3,NT"3), ATPL(3,1), STPC(NT*3,NT),SV(N) 
STPS(NT*3,NI*3), STPL(NT*3,1), L1(NOBS,1), L(NT,1), 
ICOL(NV) , IROW(NV), SPS(NT*3,NT*¥3),SPLC(NT*3, 1), ANORM(NN), 
XHAT(NP),XACNP) ,VTP(1,NT),V(NT,1),V1(NOBS),SIG(1,1), @ 
SDX(NT) , SDY(NT) , SDZ(NT) , P2(NP) , L2(NP) , LB2( NP) , AUX(2) 


NOTE: THE SAME TRANSPONDERS MUST APPEAR IN EVERY PEGASUS 
POSITION. A RECORD IN WHICH ONE (OR MORE) DROP OUT 


@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
@ 
IS NOT ALLOWED. @ 
@ 


@AAEECAEAERAAAEEERAAECEAARAEAAAARCARCEEEAEEAREAEAAEEEAREEECEGAEEGEEE 


DIMENSION XO( 72),P(4,4),BIGP(4, 80),TRANS(4) ,X(4),Y(4),2(4), 

1 A(4,3),8(4,12),VALUE(480 ),ATP(4,3) ,ATPA(3,3),ATPS(3,12),SV(20), 
2 ATPL(3,1), STP(12,4), STPS(12,12), STPL( 12,1), ANORM( 2628), 
3 XHAT( 72),XAC 72), VTP(1,4),V(4,1),V1( 80),SIG(1,1), SDX(4), 
4 SDY(4), SDZ(4), P2( 72), AUX(2) 

DOUBLE PRECISION LB( 80,1), L1( 80,1), L(4,1), L2( 72), LB2( 72) 
CHARACTER*12 PROJ 

INTEGER*2 ICOL(480 ), IROW(480 ) 


COMMON SPS(12,12), SPL(12,1) 


AEEAEEEEEECEEAEEEEEEAAEECEEEEEEEECEEECEECEEEECEEAEEECEECCEEEEEEEEE 


@ @ 
@ READ IN THE FiXED DATA AND ANY ASSOCIATED VARIANCES @ 


gFS, 


C)Glrca 2° G2) 02 


CG) OC) GaG@ee) cB 


Cie) 


@ UNITS 5 AND 6 READ AND WRITE FROM AND TO THE SCREEN G 
@ FOR INTERACTIVE PROCESSING. UNITS 1 AND 2 READ AND @ 
@ WRITE FROM AND TO DISK FILES. @ 
@ @ 
CECLECAECECAECAECEEEEEEEEEEECEECCEEEEEEECEEECEEEEEECCCECEEEEEEEEG 


WRITE(6,*) 'PROJECT OR RUN DESCRIPTION (TO 60 CHARACTERS) ' 
READ( 1,10) PROJ 
10 FORMAT(A12) 
WRITE(6,*) ‘NUMBER OF TRANSPONDERS (THEIR COORDINATES AND STD. ' 
1 ‘DEVIATIONS TO BE READ FROM THE DISK FILE) ' 
READ(1,*) NT 
READ( 1%) C81); YG), 201) ee 
READ(1,*) CSDX(1),SBYCL)..SDZCD ie 
WRITE(6,*) ‘THE GLOBAL STD. DEVIATION TO BE GIVEN TO THE ', 
1 'Z PEGASUS COORDINATES’ 
READ(1,*) SDZP 
WRITE(6,*) ‘CONVERGENCE LIMIT FOR ADJUSTMENT AND # OF ITERATIONS’ 
READ(1,*) TOL, NITER 
WRITE(2,15) PROJ,NT,(X(1),SDX(1),YC1),SDY(1),Z€1),SDZ(1), 1 =i 
15 FORMAT('1' ,///,5X,70@°* ),// /UOX,A60 9/7/ =oxe0e = le 
1 'NUMBER OF TRANSPONDERS =' ,12,//,'TRANSPONDER COORDINATES AND ‘|, 
2 'THEIR STANDARD DEVIATIONS’ ,//,4(10X,3(F10. 2,2X,F6.2),/),///) 
WRITE(2,16) SDZP 
16 FORMAT(' ',5X,'THE GLOBAL STANDARD DEVIATION FOR THE PEGASUS ‘, 
1 'Z COORDINATES =' ,F7. 2) 
WRITEC2 417). EOL enone: 
17 FORMAT(' ',//,5X,'THE CONVERGENCE LIMIT FOR THE ADJUSTMENT =' , 
1 F7.3,//,6X,'THE NUMBER OF ITERATIONS ALLOWED =', 13,///, 
2 50X, INPUT DATA’ ,//>5X%,, PRESS 0.) eee2e TIME DELAY 1 Bee 
3 'TIME DELAY 2' ,2X,‘° TIME DELAY 3' ,2X, TIME DELAY 4°, 5X, S@UNaae 
4° VELOCITY’ , 2X, APPROX. PEGASUS POSITIONS(X, Yaz) ,7/) 


-aeiebiebbbeiiiabii Dea ceaa a ag ESSEC EES 
@ READ IN THE PEGASUS DATA, POSITION BY POSITION @ 


@ @ 
GEEECEEEEEEEEEECEECEEEEEECEEEEECEEEECEEEEEECEEEEEEEEEEEEEECECEEEE 


FIRST ZERO OUT THE P MATRIX 


b 


DO 20 1 = 4 NE 

BOe20 Mies lena 
PCr) e=" 0. ODO 

20 CONTINUE 
N= 0 

25 READ(1,*, END = 45) PRESS, ZP,(TRANS(I), I=1,NT),SV(N+1) 
NNT = N*NT 
COs One —alena 
LB(NNT + I,1) = TRANS(I) 

30 CONTINUE 


CALCULAT® THE APPROXIMATE POSITION OF PEGASUS 
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Gre) OO @ 2 1 Gr 


C2 G2 ©) 


O10) 2 Gren) 2 OG are 


CALL APPRO(TRANS,X,Y,Z,SV(N+1) ,NT, AUX, XP, YP, ZP) 


XO(N*3 + 1) = XP 
XO(N*3 + 2) = YP 
XO(N*3 + 3) = Z2P 
WRITE( 2,32) PRESS,2ZP,(TRANS(1I),I=1,NT),SV(N+1),XP,YP,ZP 


32 FORMAT(' ',3X,F7.2,2X.F7. 2,5X,F7.4,3(6X,F7.4) , 13X,F8.3, 7X, 
Meee 7. ))) 


DO ERROR PROPAGATION AND FORM THE P MATRIX FOR THE FIRST PEGASUS 
POSITION. STORE THIS IN MATRIX BIGP. THE ERROR PROPAGATION 
SUBROUTINE MUST BE INSERTED HERE AND THE NEXT FEW LINES OF CODE 
MODIFIED ACCORDINGLY. 


EE&EEEEEEEEEEEE6EEEE6666hE666666666666666565 
FOR INITIAL PROGRAM TESTING, INSERT DIAGONAL VALUES ONLY 


DOe 35 le — NT 
Pol 1 ))= 1 000/70. 0005D0** 2 
35 CONTINUE 
&&E6666Eh566E656666566666666666666666666666& 


DO 40 I seedy 
DO 40 J ene 
BIGE Gt NNT tJ) = PCI ,J) 
40 CONTINUE 
N=N+1 
GO TO 25 
45 CONTINUE 


CALCULATE PROGRAM CONSTANTS 


NMAX = N 
NP = N*3 + (NT*3) 
NOBS = N*NT 


NV = (6*NT) * N 
NN = NP*(NP+1)/2 
NT3= NT*3 
ITER = 0 

WRITE( 2,46) NMAX,NOBS,NP 


46 FORMAT('1',///,5X,'ACTUAL ## OF PEGASUS POSITIONS =',15,//, 

i Sn TODAL 7 GF OBS. ON PEGASUS = ,15,//, 

2 5X, TOTAL # OF PARAMETERS =',15) 
E@EECEEECEEEEEEEECEEEEEEEEECEEEEEEEEEEEEEECCCEECEEECEEGEAEECEAEEECGGE 


SET UP THE MATRICES NEEDED FOR THE ADJUSTMENT, BEGINNING 
WITH THOSE ASSOCIATED WITH THE "OBSERVED" PARAMETERS. 
HOS Or THESE NEED TO BE ZEROEB OUT AT THIS POINT. 


THEN SET UP THE NORMAL MATRIX BLOCK BY BLOCK, EACH BLOCK 
CORRESPONDING TO A NEW PEGASUS POSITION 


CEEEEEEEACEEECEEECEEEECEECEEEEEEEEEEEEECECEEECEEEEEEEEEEEEEEEEEEG 


DAADADADAMMAM 
WAOAODADDADM|AM 
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Slee: 


ep) ap a 


Gi G CC) C2 ct) Ca 


ZERO OUT MATRICES 


DO 50 I = 1,NP 
LB2(1) = 0. ODO 


E2( 1) = 0. 0D6 

L2(1) = 0. 0bD0 

XHATC 1) = 0. ODO 
50 CONTINUE 


INSERT DHE NON-ZERO SEEEMEN IS 


J = N*3 
DO 55 I = 3,J,3 
P2(1) = 1. ODO/SDZP**2 
LB2(1)= XO(TI) 

55 CONTINUE 


J =N*3 + 1 

K=1 

DO 58 I = J,NP,3 

XOC1) = X(K) 

P2(1) = 1. ODO/SDX(K)**2 
LB2(1) = XO(T) 

XO(I+1) = Y(K) 

P2(I+1) = 1. ODO/SDY(K)**2 
LB2(1+1)= XO(I+1) 

XO(1+2) = 2(K) 

P2(1+2) = 1. ODO/SDZ(K)**2 
LB2(1+2)= XO(1+2) 
K=K+1 


58 CONTINUE 
WRITE( 2 ,800)(LB(I,1),1=1,NOBS) 

800 FORMAT( ‘1° ,///,50X, LB MATRIX ,//,8(C1X, 10(F10. 7,2X)/), 1X, ZF poe 
WRITE( 2 ,801)(XO(1) ,1=1,NP) 

801 FORMAT( 0° ,///,50X,'XO VECTOR’ ,//,8(C1X, 10CP10. 4, 2x)7 tx RiGee 
WRITE(2,802)(P2(1) ,I=1,NP) 

802 FORMAT('0O',///,50X,'P2 VECTOR’ ,//,8(1X,10(F10. 4,2x)/ 1k. Blom 
WRITE( 2 ,803)(LB2(1),1=1,NP) 

803 FORMAT( '0',///,50X,'LB2 VECTOR ,//,8C(1X, 10(F10, 4, 2X)/) 1%. F me 


Be re aS Oa 


@ 
@ THE ITERATIVE PROCESS BEGINS FROM HERE. BEGIN BY ZEROING ea 
@ OUT THE NORMAL MATRIX AND THE TWO COMMON BLOCK MATRICES @ 


@ @ 
AEEEEEEEEEECEEEEEEEEEEEEEEEEEEEEEEEEEEEECEEEEEEEEEECEEEEEEEEEEEEC 


60 DO 62 I = 1,NT3 
SEL@a = 07000 
DO 62 J = 1,NT3 
SPS(I,J) = 0. 0DO 
62 CONTINUE 
DO 65 I = 1,NN~ 
ANORM(1) = 0. ODO 
65 CONTINUE 


SZ 


Care) C2 C2 


Gi G2 ©) 


Ci-G)C 


kK = 1 


J Nei eat 
DOToo ie ==) Nr 3 
X(K) = XOCT) 
YCK) = XOC1I+1) 
ZCK) = XOCI+2) 
K = Ktl 
66 CONTINUE 


BEGIN THE FORMATION OF THE A AND Ll MATRICES 


68 


70 


le fSVCK)) 


DO 100 K = 1,NMAX 

[= OS ash 

DO 68 I = 1,NT 

TRANS(I) = LB(JJ + 1,1) 
DO 68 J = 1,NT 

P(I,J) = BIGP(I,JJ + J) 
CONTINUE 

KK = (K-1)*3 

re OC KK er) 

YP = XO(KK + 2) 

ZP = XO(KK + 3) 


nou ue ul 


DO 70 I = 1,NT 

Li(JJ + 1,1) = (DSQRT((XP-X(1))**2 + (YP-Y(I))**2 + (ZP-Z(1))**2) 
- TRANS(I) 

ene) c= ea@idects 1.1) 


CONTINUE 


ZERO OUT THE A AND S SUBMATRICES 


75 


80 


DO 80 I = 1,NT 
DO ee) = 136 
ACI,J) = 0. 0D0 
CONTINUE 


DO 80 J = 1,NT3 
S(I,J) = 0. ODO 
CONTINUE 


CALL SETUP(X,Y,Z,XP,YP,ZP,SV(K) ,K,NT,NV,NT3,NMAX,A,S,ICOL, IROW, 
1 VALUE) 


CALL NORM(A,S,P,L,K,NT,NT3,NN,NP,NMAX,ATP,ATPA,ATPS, ATPL, STP, 
1 STPS,STPL,ANORM, XHAT) 


100 CONTINUE 


WRITE(2 ,804)(L1¢1,1) ,1=1,NOBS) 


SOGerORMAMC sihey/ (250%, Lie MATRIX: .//58C1X,10(F 10. 7,2X)/),1X,2P10. 7) 


PLACE SPS AND SPL INTO THE NORMAL MATRIX AND THE COLUMN VECTOR 
RE SOPhe ily rao 


DO 110 I = 1,NT3 
NCGin = CNP + 1)ce "1 
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Ci C2 C2 Ga CO CuCa Co to 


mC) 


Se) BP) 


Ii = NCOL*(NCOL + 1)/2 + 1 

WO we a = 1 Nee 

ee 

ANORM( 12) = SPS(NT3-J+1,NT3-I+1) 
110 CONTINUE 

13 = NMAX*3 + 1 

DO 120 I = 13,NP 

XHAT( 1) = SPLO1- soni) 
120 CONTINUE 


SET UP THE L2 MATRIX AND ADJUST THE DIAGONAL ELEMENTS OF THE 
NORMAL MATRIX FOR THE INFLUENCE OF P2 


CALL P2L2(P2,LB2,X0,ANORM,L2, XHAT ,NP,NN) 


dab ESOL G02 SOS ESE LEEEE SSE SECT 
@ 
@ SOLVE THE NORMAL EQUATIONS AND COMPUTE THE RESIDUALS @ 


@ @ 
CEEEEECEEEEEEEEEEEEEEECEEEEECEECECECEEEECEEEEEEEEEEEEEEEEEEECE 


CALL A USER SUPPLIED ROUTINE TO SOLVE THE NORMAL EQUATIONS. 
(SEE APPENDIX A, SECTION D.) 


DO 135 I = 1,NP 

XHATCI) = -XHAT(I) 

XAC1) = xOCr ee Ana) 

XO ec) 
135 CONTINUE © 

ITER = ITER + 1 

WRITE( 2,140) 
140 FORMAT('0' ,///,5X,' ADJUSTED PEGASUS POSITIONS (WITH THE '‘, 

1 'CORRECTIONS TO THE APPROX. POSITIONS IN PARENTHESES)',//, 18X, 

2: Xe ABR Yoel eae) 

J = NP - NT3 

DOniees iene 

WRITE( 2,142) XACI),XHATCI),XACI+1),XHAT(I4+1) ,XA(I+2), 

1 XHAT(I+2) 
142 FORMATC’ °,10X,F9.3.2X,'(',Fés1, |) 5202% F9.3,2% 510 hoe 
145 CONTINUE 

DF Satin 

WRITE( 2, 146) 
146 FORMAT(' ',//,5X,'ADJUSTED TRANSPONDER POSITIONS (WITH', 

1 'CORRECTIONS TO THE APPROX. POSITIONS IN PARENTHESES)’ ) 

DO 148 I=J,NP,3 

WRITE( 2,142) XAC(I),XHATCI),XACI4+1) ,XHATCI+1) ,XAC1+2), 

1 XHAT(I+2) 
148 CONTINUE 


TEST THE PARAMETER CORRECTIONS TO SEE IF ITERATION IS REQUIRED 
IFLAG = 0 
DOGS 0 sie =1 (NE 


IFC DABS(XHAT(I)).GT. TOL) IFLAG = IFLAG + 1 
150 CONTINUE 
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o) 


Gr ca Gs CG) CarPC) Gl C3 Gra Came Co C2 Gira 


OP Se Fal Se i > LB 


Biel GH ys Sel 


COMMENT OUT THE NEXT THREE STATEMENTS DURING PROGRAM DEVELOPMENT 
TE CUPEAGZEQ: 0) GO TO 155 


TPGITER SE NITER) GO 10°60 
COMPUTE. THE RESIDUALS 


155 CALL RESID( VALUE, ICOL, IROW, BIGP,L1,XHAT,NMAX,NT,NV,NP,NOBS, 
Py EP G51, 5160) 


WRITE(2, 160) 
160 FORMAT( 1',///,10X,' OBSERVATIONAL RESIDUALS AFTER THE ADJUSTMENT’ , 
1 //,5X,' TRANSPONDER 1 TRANSPONDER 2 TRANSPONDER 3° ,4X, 


2' TRANSPONDER 4' ) 

DO 170 I = 1,NOBS,NT 

WRITE(2,165)(V1(I+K-1), hoe eel) 
165 FORMAT( ',9X,F8.4,3(9X,F8. 4)) 
170 CONTINUE 


WRITE(2,175) SIGO 
175 FORMAT(' ',///,5X,'THE A POSTERIORI VARIANCE OF UNIT WEIGHT =', 
1 F8.4) 


*THE FOLLOWING WRITE STATEMENT WAS ADDED BY M. HASKELL TO LOOK 

"AT THE VARLANCE-COVARTANCE MATRIX FOR POSITION COORDINATES. 

*THE CALL STATEMENT FOR SUBROUTINE MAT WAS ADDED TO ACCESS ONLY 
Hee iNeOR iat LON NEEDED TO CALCULATE THE VARTANCES AND COVARIANCES 
SOP sii VELOCITIES. 


WRITE(2,*)(ANORM(II),II=1,NN) *PROVIDES (ATPA)7~! IN vECTOR FORM 
CALL MAT( ANORM ,NMAX) 


Pi eer races ccc coc ec cocec cc roerr or cecee 


@ COMPUTE VARIOUS STATISTICAL QUANTITIES INCLUDING THE @ 
@ VELOCITY OF PEGASUS. @ 
C 


@ 
CEEECECEEEEEAEEEEEEEEEEEEEEEEEEEEEECEEEECEEEEEEEEEEEEEEEEEEECEG 


SLOP 
END 


Pith eee Coe ce cee eee 
@ 
@ SUBROUTINES USED IN THE PROGRAM @ 


@ @ 
EEEEEEECEECEEEEEEECEECECEEEEEEEEEEECEEECEEEECEEEEE 


SURROUMING A PPROGEIRANS  4.1,2,07,NT,AUX,XPlYP,ZP) 
IMPLICIT REAL*8(A-H,O-Z) 


THIS SUBROUTINE COMPUTES THE APPROXIMATE COORDINATES FOR PEGASUS 


INPUT: TRANS - VECTOR OF TRANSPONDER TIME DELAYS FOR THE DESIRED 
PEGASUS POSITION. 
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C01) Gare 


O 


GGG) ©) GG OC) C) oO.) (a 


X,Y,Z - ARRAYS CONTAINING THE COORDINATES OF THE TRANSPONDERS 
SV - VELOCITY OF SOUND IN WATER 
NT - NUMBER OF TRANSPONDERS 

OUTPUT: XP,YP,2P - COORDINATE SS@EernG.ouls 


DIMENSION TRANS(NT), XCNT),YCNT),ZCNT) ,AUX(2) 


DO 10 I = 1,NT 
TRANS(I) = TRANS(I)*SV 
10 CONTINUE 
DO 20 I = 1,2 
AUX(1) = (TRANS(I+1)**2 - (ZP-2(14+1))**2 - X(I+1)**2 - YC I+1)**2) 
1 - (TRANS(I)**2 - (ZP-ZC1))%*2 - X(1)%**2 - Y(1)**2) 
20 CONTINUE 
YP = 0. S5DO*(AUX(1)/(X(1)-X(2)) - AUX(2)/CKC2)-xX(3)))7 (De 
1 /(XC1)-XC2))e= CY(2)-YOs 7 eG eae 
XP = 0. 5D0*( AUXGI) AGG) -1 G2 = AUX(2)/C¥(2)- ¥(3))) /CCS( 1) as 
TY 7 ONG ve2)) - CXC2) - XC GG) 
DO 30 I = 1,NT 
TRANS(I) = TRANS(1I)/SV 
30 CONTINUE 
RETURN 
END 


SUBROUTINE SETUP(X,Y,Z2,4P,YP,2P,SV,N,NIT,NV,NI3,NMAX,A,S, [COLI RORs 
1 VALUE) 
IMPLICIT REAL*8( A-H,O-Z) 


THIS MATRIX SETS UP THE A AND S MATRICES NEEDED FOR THE BLOCK 
BY BLOCK FORMATION OF THE NORMAL MATRIX. 


INPUT: X(1I), I = 1,NT X COORDINATES FOR THE TRANSPONDERS 
YC) ee" NT YY “ 
ZGL) I = ies NT 7 tt 11 ti tl 
XP,YP,ZP, APPROXIMATE COORDINATES FOR PEGASUS. 
SV = VELOCITY OF SOUND THROUGH WATER. 
N = NTH POSITION OF PEGASUS. 


NT, NV, NT3 SAME AS IN THE MAIN PROGRAM. 
NMAX = ACTUAL # OF PEGASUS POSITIONS IN THIS DROP 
OUTPUT: A = THE PORTION OF THE A MATRIX CORRESPONDING TO THE 
NTH PEGASUS POSITION. 

S = AS ABOVE, BUT THE OUTER BAND OF THE A MATRIX. 

VALUE = THE NON ZERO # IN THE A MATRIX CORRESPONDING 
TO THE COLUMN # HELD IN ICOL AND THE ROW # HELD 
IN IROW. 


DIMENSION X(NT), YCNT), Z(NT), ACNT,3), SCNT,NT3), VALUE(NV) 
INTEGER*2 ICOL(NV), IROW(NV) 


NROW = (N-1)*NT 
NCOL = (N-1)*3 
NCOL1 = NMAX*3 
ICOUNT = (N-1)*6*NT 
DO 20 I = 1,NT ° 
: = 1 
= ([- 1)*3 + ]j 
coe DSQRTCCXP: = XC 1) #2 Cee) 2 CZ eas 
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NOOO Oa ea ae ae 


ACJ) = (XPe= X01))/(DIST*SV) 
ICOUNT = ICOUNT + 1 
VALUE(ICOUNT) = A(I,J) 

TROW( ICOUNT) = NROW + I 
ICOL(ICOUNT) = NCOL + J 
A(I,J+1) = (YP - Y(I))/(DIST*SV) 
ICOUNT = ICOUNT + 1 
VALUE(ICOUNT) = A(I,J+1) 

IROW( ICOUNT) = NROW + I 

ICOL( ICOUNT) = NCOL + J+1 

ACI ,J+2) = (ZP - 2(1))/(DIST*SV) 
ICOUNT = ICOUNT + 1 
VALUE(ICOUNT) = A(I,J+2) 
IROW(ICOUNT) = NROW + I 

REOL( COUNT) w= NCOL 4) J+2 
S(1l-K) = -AC Ped) 

ICOUNT = ICOUNT + 1 
VALUE(ICOUNT) = S(I,K) 
IROW(ICOUNT) = NROW + I 
ICOL(ICOUNT) = NCOL1 + K 

S(T Kil) = =AC1,J+1) 

ICOUNT = ICOUNT + 1 
VALUE(ICOUNT) = S(I,K+1) 

IROW( ICOUNT) = NROW + I 
ICOL(ICOUNT) = NCOL1 + K+l1 
S(I,K+2) = -A(I,J+2) 

ICOUNT = ICOUNT + 1 

VALUE( ICOUNT) = S(I,K+2) 


IROW(ICOUNT) = NROW + I 

ICOL(ICOUNT) = NCOL1 + K+2 
20 CONTINUE 

WRITE(2,100) 


100 FORMAT('1',///,15X,'A MATRIX AND THEN THE S MATRIX, BLOCK BY', 
lee SB LOCK " ) 
WRITE(2,101)((ACI,J) ,J=1,3) ,1=1,4) 
101 FORMAT(' ', T2,4(2X,3(F10.8,3X),/)) 
WRITE( 2 ,102)((S(1,J),J=1,12) ,1=1,4) 
102 FORMAT(' ', T2,4(12(F10.7,1X),/)) 
RETURN 
END 


SUBROUTINE NORM(A,S,P,L,N,NT,NT3,NN,NP,NMAX,ATP,ATPA,ATPS,ATPL, 
1 STP ,STPS ,STPL,ANORM, XHAT) 
IMPLICIT REAL*8(A-H,0-Z) 


THIS SUBROUTINE FORMS THE NORMAL EQUATIONS, BLOCK BY BLOCK 

INPUT: A MATRIX FROM SUBROUTINE SETUP 
S) 
P THE WEIGHT MATRIX CORRESPONDING TO THIS BLOCK 
ete eA Trix 
Ne NE, ONTS, AS. FOR SUBROUTINE SETUP 
NN, NP, NMAX , AS FOR MAIN PROGRAM 

OUTPUT: A NUMBER OF AUXILLARY MATRICES USED IN THE MATRIX 
NUDE EICATIONS, le AIP ATPA, ALPS, AITPL,STP, SIPS, 
AND STPL. 
MATRICES SPS AND SPL WHICH COME IN VIA THE COMMON 
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ap TE @ a SP TS 


Lee pe 


De Se Ne Ie 


2) Cane 


BLOCK ARE UPDATED. 
ANORM - THE NORMAL MATRIX (CIN ARRAY FORM) 
XHAT - THE ATPL VECTOR 


DIMENSION ACNT,3), SCNT,NT3),PCNT, NI) ADECND 33). Ate eGs on. 


1 ATPS(3,NT3), ATPL(3,1), STP(NT3,NT), STPS(NT3,NT3) ,STPL(NT3,1), 
2 ANORM(NN) ,XHAT(NP) 


DOUBLE PRECISION L(NT, 1) 
COMMON SPS(12,12)5 SPLU UZ) 


PERFORM THE MATRIX MULTIPLICATIONS 


100 


110 


CALL ATB(A,P,ATP,NT, 3,NT) 

CALL ABCATE. AA DPAR 2 Nias) 

CALL AB(ATP,S,ATPS,3,NT,NT3) 

CALL ABCATE, L ATPL 2 Nt 

CALL ATB(S,P,STP,NT,NT3,NT) 

CALL AB(STP,S,STPS,NT3,NT,NT3) 

CALL) ABCSTIP,L,S1PL NES Ni 

WRITE( 2, 100)(CATPA(I, J) ,J=1,3) ,ATPL(I,1),1=1, 3) 
FORMAT( ',/,15X, '‘ATPA BLOCK’ ,10X,'ATPL VECTOR’ ,//, 


l 302%, 3(D 10,352) ex, On ae 


WRITE( 2 
FORMAT( 


110)((STPS(I,J),J=1,12),1=1,12) 


tr 7 -15X, 'STPS BLOCK',//, 


V2 121 2C Diese) 


PLACE ATPA AND ATPL IN THEIR APPROPRIATE POSITION IN EITHER THE 
NORMAL MATRIX OR THE COLUMN VECTOR 


10 


20 


NGOL — Guede cueee 

11 = NCOL*(NCOL + 1)/2 
ANORM(I1) = ATPA(1,1) 
Ill = Il + NCOL 


-ANORM(T11) = ATPAC1,2) 


ANORM(I11 + 1) = ATPA(2,2) 
Tid. = Tiere | Noor 
ANORM(1111) = ATPA(1,3) 
ANORM(I111 + 1) = ATPA(2,3) 
ANORM(I111 + 2) = ATPA(3,3) 


DO 10 I = 1,NT3 

INT = (NMAX*3) + I 

J = INT(INT - 1)/2 + 1 + (N-1)*3 
ANORM(J) = ATPS(1,T) 

ANORM( J+1) = ATPS(2,1) 

ANORM( J+2) = ATPS(3,1) 

CONTINUE 


DO 20 I = 1,3 
XHAT(NCOL+I-1) =ATPL(I,1) 
CONTINUE 


UPDATE SPS AND SPL 


DO 30 I = 1,NT3 
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7G) C2) © 2 Cate G2) Gy CD 


Cl GO) Cy oy GG 


Sle Si Lael) + STPL( 1,1) 

Dors0 J°="15Nr3 

Soe) = Sr SG) jeer SIPS(1,J) 
30 CONTINUE 

RETURN 

END 


SUBROUTINE ATBCA,B,R,L,M,N) 
Pi Eter REAL SCA 0-Z } 


THIS SUBROUTINE COMPUTES THE MATRIX PRODUCT A’B 
CNT OT: MATRIX A (L X M) 

MATRIX B (L X N) 
OUTPUT MATRIX R (M X N) 


DIMENSION A(L,M),B(L,N),R(M,N) 
DO 10 I = 1,M 
DO 10 J = 1,N 
R(I,J) = 0.0D0 
DO 5 K = 1,1 
RG = RCI) wath ol ye ECh.) 
5 CONTINUE 
10 CONTINUE 
RETURN 
END 


SUBROUTINE AB(CA,B,R,L,M,N) 
IMPLICIT REAL*8(A-H,O0-Z) 


FORM THE MATRIX PRODUCT R = AB. 
THE MATRICES A AND B ARE RETURNED UNCHANGED 


DIMENSION ACL,M),BCM,N),RCL,N) 


DO5 I = 1,L 

DO 5 J = 1,N 

R(I,J) = 0.0D0 

DO 5 K = 1,M 

R(I,J) = R(I,J) + ACI,K)*B(K,J) 
5 CONTINUE 

RETURN 

END 


SUBROUTINE RESID( VALUE, ICOL, IROW,BIGP,L1,XHAT,NMAX,NT,NV,NP,NOBS, 
ie PW ology oGO)) 
IMPLICIT REAL*8(A-H,0-Z) 


COMPUTES THE RESIDUALS ON ALL THE TRANSPONDER TIME DELAY 
OBSERVATIONS AND THE A POSTERIORI VARIANCE OF UNIT WEIGHT. 
INPUT: VALUE, ICOL,IROW,- MATRICES DERIVED IN SUBROUTINE SETUP 
BIGP - THE FULL WEIGHT MATRIX FOR ALL OBSERVATIONS 
lest - THE VECTOR OF "COMPUTED-OBSERVED' OBSERVATIONS 
XHAT - THE LEAST SQUARES SOLUTION VECTOR 
NMAX,NT,NV,NP,NOBS - AS IN THE MAIN PROGRAM 
AUXILLARY MATRICES : P, VTP, V,SIG 
OUTPUT: V1 - THE VECTOR OF RESIDUALS 
SIGO - THE A POSTERIORI VARIANCE OF UNIT WEIGHT 
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10 


15 


16 


18 


20 


30 


40 


DIMENSION VALUE(NV) , XHATCNP) ,V1(NOBS) , BIGP(NT,NOBS) ,PCNT,NT), 


1 VTP(1,NT),SIG(1,1),VCNT, 1) 


DOUBLE PRECISION L1(NOBS) 
INTEGER*2 ICOLCNV), IROWCNV) 


DO 10 I = 1,NOBS 
V1iCI) = 0.0D0 
CONTINUE 
Ic = 1] 
DO 20 I = 1,NV 
IF( IROW(1).NE.IC) GO TO 18 


V1(IC) = V1CIC) + VALUE(1)*XHAT( ICOL(I1)) 
GO TO 20 
V1(IC) = V1CIC) + LiI(1¢C) 
ic eel 
IF(IC. LE.NOBS) GO TO 16 
CONTINUE 
V1(NOBS) = V1(NOBS) + L1(NOBS) 
SiGe 07 000 
Ic = 0 


DO 40 I = 1,NMAX 
DO 30 J = 1,NT 
V(J,1) = VICIC*NT + J) 
DO 30 K = 1,NT 
P(J,K) = BIGP(J,IC*NT + K) 
CONTINUE 
CALL ATB(V,P,VTP,NT,1,NT) 
CALL AB(VIP,V,SIG,1,NT,1) 
SIGO = SIGOm one) 
iG. = [Cm 
CONTINUE 
SIGO = SIGO/FLOAT(NOBS-NMAX*3 ) 
RETURN 
END 


SUBROUTINE P2L2(P2,LB2,X0,ANORM, L2,XHAT,NP,NN) 
IMPLICIT REAL*8(A-H,0-Z) 


THIS SUBROUTINE COMPUTES THE L2 MATRIX AND ADDS IT TO THE COLUMN 
VECTOR. IN ADDITION IT INCREMENTS THE DIAGONAL ELEMENTS OF THE 
NORMAL MATRIX TO INCLUDE THE INFLUENCE OF P2 


DIMENSION P2(NP), XO(NP),ANORM(NN) , XHATCNP, 1) 
DOUBLE PRECISION L2(NP), LB2(NP) 


DO 10 I = 1,NP 

WAI) = NOS) = JAC 

INT = I*(I+1)/2 

ANORM(INT) = ANORM(INT) + P2(1) 
P2L, = seal) eee) 

XHAT(I,1) = XHAT(1,1) + P2L 
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TO © 


HO CONTINUE 
RETURN 
END 


SUBROUTINE MAT(C,NMAX) 

*NOT PART OF J HANNAH'S ORIGINAL PROGRAM. ADDED BY M. HASKELL. 
THIS PROGRAM READS APPROPRIATE VALUES FROM THE VARIANCE-COVARIANCE 
MATRIX FILE (ANORM) AND PLACES THEM IN THE MATRICES NEEDED TO COM- 
PUTE VARIANCES OF THE CURRENT VELOCITIES. 

IMPLICIT REAL*8 (A-H,O-Z) 

DOUBLE PRECISION C(1) 

DINENSION C(444153) 

I(K) = K*(K+1)/2 

= 1 

ROW = K 

WRITE (2,*)C(I(K)), CCI(K+1)-1), CCIC(K+2)-2) 

WRITE (2,*)C(I(K+1)-1), CCI(Kt1)), CC IC K+2)-1) 

WRITE (2,*)CCI(K+2)-2), CCI(K+2)-1), CCIC(K+2)) 

DO 15 J = 1,NMAX+1 
WRITE (2,*)CC1CKH+3)-3), GCICK+4)-4), CCI(K+5)-5) 
WRITE (2,%)CCI(K+3)-2), CCICK+4)-3), CCI(K+5)-4) 
Wiitee( 2S yOGlGhta jel). —©C1L( hrs )=2), CC 1(K+5 )-3) 


Rone ae OVO CLGIEES ne COLCKES i- 1.) CC LC K+5)-2) 
WRITE (2,*)CCI(K4+4)-1), CCI(K+4)), CCICK45)=1) 
WRITE (2,*)CC(I(K+5)-2), CCI(K+5)-1), CCICK+S)) 
K = K+3 
15 CONTINUE 


WRITE (2,*)CC(I(K+3)-6), CCI(K+4)-7), CCIC(K+5)-8) 
WRITE (2,*)C(I(K+3)-5), CCI(K+4)-6), C(I(K+5)-7) 
WRI EO eCiCKt3)a4) , CCICK+t4)-S) CCICKtTS) -6) 


WRITE (2,%*)C(I(K+3)-3), CCI(K+4)-4), CCICK+5)-5) 
WRITE (2,*)CC(I(K+3)-2), CCICK+4)-3), CCICK+5)-4) 
WRITE (2,*)C(ICK+3)-1), CCI(K+4)-2), CC I(K+5)-3) 


WRITE (2,*)C(I(K+3)), CCI(K+4)-1), CCI(K+5)-2) 
WRITE (2,*)C(I(K+4)-1), CCI(K+4)), CCI(K+5)-1) 
WRITE (2,*)C(I(K+5)-2), CCICK+5)-1), CCICK+5)) 


Sos | 
WRITE (2,*)C(I(K+3)-9), CCI(K+4)-10), CCI(K+5)-11) 
Hane ue e@nGkt S26), 9 CC 10K 1-9), CCI(K+5)-10) 
WRITE (2,*)C(I(K+3)-7), CCICK+4)-8), CCICK+5)-9) 


WRITE (2,*)CC(I(K+3)-6), CCI(K+4)-7), CCI(K+5)-8) 
WRITE (2,*)C(I(K+3)-5), CC 1(K44)-6), CCICK+5)-7) 
WRITE (2,*)C(I(K+3)-4), CCI(K+4)-5), CCIC(K+5)-6) 


Wie ee Oe llicKt j= 3 2) CCICKTS)=4). ClICKTS)=5) 


WRITE (2,*)CCICK+3)-2), C(I(K+4)-3), CC I(K+5)-4) 
WRITE (2,*)CCICK+3)-1), CCI(K+4)-2), CCICK+5)-3) 


9.1. 


ais) I ae Tal ip iy Jel ip 1 > Ya ap PL ap 1 >) al ep a ip 


WRITE (2,*)CCICK+3)), CC TCHS) = ere) 
WRITE (€2,*)CCICK+4)-1), COMGR+4 )) ye COni = oo 
WRITE €2,*)CCICK+5)-2))) CCICK+S )-1 ) eC Giese 


K = 36 

DORs] 
WRITE 
WRITE 
WRITE 


WRINE 
WRITE 
WRITE 


K = Kt+36 
CONTINUE 


RETURN 
END 


= 1,6 

(2, *)CCICK43)-3), CCIC Koja CClee cess 
(2,*)CCICK+3)=2), CCICKESe) CG a) 
(2,*)C(1(K43)-1), CC1CK+a9- 2) Chas) ee) 


(2,*)CCICKt3)), CONCK+G) - 1) ee) 


(2,*)CCI(K+4)-1), CCICK+4)), CCICK+5)i=1) 
(2,*)CCICK+5)-2), CCICK+5)-1), CCICK+5)) 


22 


APPENDIX B 


CURRENT VELOCITY VARIANCE-COVARIANCE PROGRAM 


This program computes the variance-covariance matrix of 
the current velocities by a procedure described in Chapter VI, 
Section D. It uses data from the Pegasus position variance- 
covariance matrix produced by the Fortran program Pegasus as 


described in Appendix A. 
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THIS PROGRAN USES DATA FROM THE VARIANCE-COVARIANCE MATRIX PRODUCED 
BY THE PEGASUS FORTRAN PROGRAM. IT COMPUTES THE VARIANCE-COVARI - 
ANCE MATRIX OF THE CURRENT VELOCITIES. 


fa GC 1G. 


IMPLICIT REAL*8(A-H,0-Z) 
REAL SP103,3),.SP12(3.,3)5 SB26323)39 G59) OUR en 


N=1 
DO 100 I = 1,3 
READ(1,*)(SP1(1,J),J = 1,3) 
C WRITE(2,120)(SP1(I,J),J = 1,3) 
120 FORMAT(' ',3X,F10.4,3X,F10. 4,3X,F10. 4) 
100 CONTINUE 


10 DO 200 I = 1,3 
READ(1,*,END=45)(SP12(1,J),J = 1,3) 
C WRITE(2,120)(SP12(I,J),J = 1,3) 
200 CONTINUE 
DO 300 I = 1,3 
READ(1,*)(SP2(1,J),J 
¢ WRITE(2,120)(SP2(I,J), 
300 CONTINUE 


= 1,3) 
J = 1,3) 
E 
eG WRITE(2,15)N 
15 FORMAT( ',/,5X, | VELOCITY VARIANCE-COVARIANCE MATRIX FOR ', 
1'POSITION ~ ,13,/) 
DO 00s s— 1-3 
DO 350 J = 1,3 
V(1,J)=CCSP1C1,J)+SP201 J) -SP1201 )-SP1Z2e 7 GliGees) 
350 CONTINUE 
C WRITE(2, 120)0VGlJ).4 = les) 
400 CONTINUE 
SU = SQRTC(V(1,1)) 
SV = SQRT(V(2,2)) 
SW = SQRT(V(3,3)) 
CC WRITE(2,130)SU, SV, SW 
CC130 FORNAT( °,/,3X, SIGMA U = ,F10.4, 3X, SIGHASY = (rillgnce one 
CC 1 ‘SIGMA W =',F10.4/) 
WRITE( 2, 600)N,SU,SV,SW 
600 FORMATC ~,10X, POS 13,3, 8724530 7 oe) 
DO™%50 = 3 


. DO 460 J = 1,3 
SP1(1,J) = SP2(1,J) 
460 CONTINUE 
C WRITE(2,140)(SP1(I,J),J=1, 3) 


C 140  FORMAT(! 


450 CONTINUE 
N=N +1 
GO TO 10 

45 CONTINUE 
STOP 
END 


‘3X, Fide dX Flo. GeoKer los) 
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